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ABSTRACT 

Over the past several decades, galaxy formation theory has met with significant suc- 
cesses. In order to test current theories thoroughly we require predictions for as yet 
unprobed regimes. To this end, we describe a new implementation of the Galform 
semi-analytic model of galaxy formation. Our motivation is the success of the model 
described by Bower et al. in explaining many aspects of galaxy formation. Despite 
this success, the Bower et al. model fails to match some observational constraints and 
certain aspects of its physical implementation are not as realistic as we would like. 
The model described in this work includes substantially updated physics, taking into 
account developments in our understanding over the past decade, and removes certain 
limiting assumptions made by this (and most other) semi-analytic models. This allows 
it to be exploited reliably in high-redshift and low mass regimes. Furthermore, we have 
performed an exhaustive search of model parameter space to find a particular set of 
model parameters which produce results in good agreement with a wide range of ob- 
servational data (luminosity functions, galaxy sizes and dynamics, clustering, colours, 
metal content) over a wide range of redshifts. This model represents a solid basis on 
which to perform calculations of galaxy formation in as yet unprobed regimes. 

Key words: galaxies: general, galaxies: formation, galaxies: evolution, galaxies: high- 
redshift, intergalactic medium 



1 INTRODUCTION 

Understanding the physics of galaxy formation has been an 
active field of study ever since it was demonstrated that 
galaxies are stellar systems external to our own Milky Way. 
Modern galaxy formation theory grew out of early studies 
of cosmology and structure formation and is set within the 
cold dark matter cosmological model and therefore proceeds 
via a fundamentally hierarchical paradigm. Observational 
evidence and theoretical expectations indicate that galaxy 
formation is an ongoing process which has been occurring 
over the vast majority of the Universe's history. The goal 
of galaxy formation theory then is to describe how under- 
lying physical principles give rise to the complicated set of 
phenomena which galaxies encompass. 

Approaches to modelling the complex and non-linear 
processes of galaxy formation fall into two broad categories: 
direct hydrodynamical simulation and semi-analytic mod- 
elling. The division is of a somewhat fuzzy nature: semi- 
analytic models frequently make use of N-body simulation 
merger trees and calibrations from simulations, while simu- 
lations themselves are forced to include semi-analytical pre- 
scriptions for sub-resolution physics. The direct simulation 
approach has the advantage of, in principle, providing pre- 
cise solutions (in the limit of large number of particles and 



assuming that numerical artifacts are kept under control), 
but require substantial investments of computing resources 
and are, at present (and for the foreseeable future), more 
fundamentally limited by our incomplete understanding of 
the various sub-resolution physical processes incorporated 
into them. The semi-analytical approach is less precise, but 
allows for rapid exploration of a wide range of galaxy prop- 
erties for large, statistically useful samples. A primary goal 
of the semi-analytic approach is to develop insights into 
the process of galaxy formation that are comprehensible in 
terms of fundamental physical processes or emergent phe- 
nomen43 

The problem is therefore one of complexity: can we tease 
out the underlying mechanisms that drive different aspects 
of galaxy formation and evolution from the numerous and 
complicated physical mechanisms at work. The key here is 
then "understanding". One can easily comprehend how a 
1 /r^ force works and can, by extrapolation, understand how 
this force applies to the billions of particles of dark matter 

^ A good example of an emergent phenomenon here is dynamical 
friction. Gravity (in the non-relativistic limit) is described entirely 
by 1/r^ forces and at this level makes no mention of frictional 
effects. The phenomenon of dynamical friction emerges from the 
interaction of large numbers of gravitating particles. 
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in an N-body simulation. However, it is not directly obvi- 
ous (at least not to these authors) how a force leads 
to the formation of complex filamentary structures and col- 
lapsed virialized objects. Instead, we have developed sim- 
plified analytic models (e.g. the Zel'dovich approximation, 
spherical top-hat collapse models etc.) which explain these 
phenomena in terms more accessible to the human intel- 
lect. It seems that this is what we must strive for in galaxy 
formation theory — a set of analytic models that we can com- 
prehend and which allow us to understand the physics and 
a complementary set of precision numerical tools to allow 
us to determine the quantitative outcomes of that physics 
(in order to make precision tests of our understanding). As 
such, it is our opinion that no set of numerical simulations 
of galaxy formation, no matter how precise, will directly re- 
sult in understanding. Instead, analytic methods, perhaps 
of an approximate nature, must always be developed (and, 
of course, checked against those numerical simulations) to 
allow us to understand galaxy formation. 

Modern semi-analytic models of galaxy formation be- 



gan with White & Frenk (19911, drawing on earlier work 
by |Rees fc Ostriker (1977|) and|White fc Rees (1978|). Since 



then numerous studies (Kauffmann et al. 199J Hole et al, 
1994 Baugh e t al. 1999I IBaugh e t al. 199^p omervill e fc PrT 
mack 199£ Cole et al. 200C' Benson et al. 2002E.Hatton et al. 



200SMonaco et al. 2007| ) have extended and improved this 
original framework. Current semi-analytic models have been 
used to investigate many aspects of galaxy formation includ- 
ing: 



• galaxy counts ( Kauffmann et al. 199^pevriendt fc 
IGuiderdoni^OOO] ); 

• galaxy clustering ([Diaferio et al. 199!^^auffmann et al.| 
1999aKauffmann ct al. 1999h Baugh et al. i999a, Benson ct al. 
'2000aBenson et al. 2000biWechsler et al. 2001Blaizot et al.i 



2006); 



• galaxy colours and metallicities ( 


Kauffmann & Char- 


lot 199^ppringel et al. 2001 


(janzoni et al. 200^Pont et al. 


200^|Nfagashima et al. 2005b 


1; 



• sub-mm and infrared (IR) galaxies (Guiderdoni et al. 
199^pranato et al. 200C|paugh et al. 200^|Lacey et al. 2008"f -' 

• abundance and properties of Local Group galaxies 



(Benson et al. 20021: jomerville 20021; 



• the reionization of the Universe (Devriendt et al, 
ll99aBenson et al. 200]IBomerville fc Livio 2003Benson et al. 



2006); 



• the heating of galactic disks (Benson et al. 2004); 

• the properties of Lyman-break galaxies (Governato 
|et al. 199^plaizot et al. 200:^plaizot et al. 2004] ); 

• supermassive black hole formation and active galactic 



nuclei (AGN) feedback ( 


Kauffmann & Haehnelt 200C Z!roton 


let al. 2006Bower et al. 


2006IVIalbon et al. 200'ilBomerville 


et al. 20081: pontanot et al. 2009a I; 


• damped Lyman-a 


systems (Mailer et al. 2001 Vlaller 



et al. 20031; 



• the X-ray properties of galaxy clusters ( [Bower et al.| 
200]|power et al. 2008| ); 



• chemical enrichment of the ICM and intergalactic 
medium (IGM) ( De Lucia et al. 200^ ^agashima et al. 



2005a I; 



galaxies (Kauffmann 199(: pe Lucia et al. 200( [^'ontanot et al. 
|200'i|pomerville et al. 2008a[). 



The goal of this approach is to provide a coherent frame- 
work within which the complex process of galaxy forma- 
tion can be studied. Recognizing that our understanding of 
galaxy formation is far from complete these models should 
not be thought of as attempting to provide a "final theory" 
of galaxy formation (although that, of course, remains the 
ultimate goal), but instead to provide a means by which 
new ideas and insights may be tested and by which quan- 
titative and observationally comparable predictions may be 
extracted in order to test current theories. 

In order for these goals to be met we must endeavour 
to improve the accuracy and precision of such models and 
to include all of the physics thought to be relevant to galaxy 
formation. The complementary approach of direct numeri- 
cal (N-body and/or hydrodynamic) simulation has the ad- 
vantage that it provides high precision, but is significantly 
limited by computing power, resulting in the need for in- 
clusion of semi-analytic recipes in such simulations. In any 
case, while a simulation of the entire Universe with infinite 
resolution would be impressive, the goal of the physicist is 
to understand Nature through relatively simple argument^ 

The most recent incarnation of the Galform model 



was described by Bower et al. (2006). The major innovation 



of that work was the inclusion of feedback from AGN which 
allowed it to produce a very good match to the observed lo- 
cal luminosity functions of galaxies. In particular, the |Bowerl 
|et al. (2006| ) model was designed to explain the phenomenon 
of "down sizing" . While the Bower et al. (2006 ) model turned 



out to also give a good match to several other datasets — 
including stellar mass functions at higher redshifts, the lumi- 
nosity function at z = 3 ( [Marchesini & van Dokkum 2007), 



the abundance of 5 < 2 < 6 galaxies ( McLure et al. 2009[ ), 
overall colour bimodality ( [Bower et al. 2006[ ), morphology 
(Parry et al. 2009), the global star formation rate and the 



black hole mass vs. bulge mass relation (Bower et al. 2006 1 — 
it fails in several other areas, such as the mass-metallicity 
relation for galaxies, the sizes of galactic disks (Gonzalez 



et al. 2009 1, the small-scale clustering amplitude (Kim et al 



2009[ ), the normalization and environmental dependence of 



galaxy colours ( Font et al. 2008 ) and the X-ray proper- 
ties of groups and clusters (Bower et al. 2010). Addition- 



ally, while the implementation of physics in semi-analytic 
models must always involve approximations, there are sev- 
eral aspects of the Bower et al. (2006) model which call 



out for improvement and updating. Chief amongst these is 
the cooling model — crucial to the implementation of AGN 
feedback — which retained assumptions about dark matter 
halo "formation" events which make implementing feedback 
physics difficult. Our motivation for this work is therefore 
to attempt to rectify these shortcomings of the [Bower et al.j 
(2006 ) model, by updating the physics of Galform, re- 



the formation histories and morphological evolution of 



^ For example, while it is clear from N-body simulations that the 
action of gravitational forces in a cold dark matter (CDM) 
universe lead to dark matter halos with approximately Navarro- 
Frenk-White (NFW) density profiles, there is a clear drive to 
provide simple, analytic models to demonstrate that we under- 
stand the underlying physics of these profiles (^Taylor &: Navarro[ 
200][parnes et al. 200"7^3arnes et al. 2007b"[ . 
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moving unnecessary assumptions and approximations, and 
adding in new physics that is thought to be important for 
galaxy formation but which has previously been neglected 
in Galform. In addition, we will systematically explore the 
available model parameter space to locate a model which 
agrees as well as possible with a wide range of observational 
constraints. 

In this current work, we describe the advances made in 
the Galform semi-analytic model over the past nine years. 
Our goal is to present a comprehensive model for galaxy 
formation that agrees as well as possible with current ex- 
perimental constraints. In future papers we will utilize this 
model to explore and explain features of the galaxy popula- 
tion through cosmic history. 

The remainder of this paper is structured as follows. In 
i 2 we describe the details of our revised Galform model. In 
[ 3 we describe how we select a suitable set of model param- 
eters. In 34]we present some basic results from our model, 
while in 95] we explore the effects of certain physical pro- 
cesses on the properties of model galaxies. Finally, in !|6]we 
discuss their implications and in ^we give our conclusions. 
Readers less interested in the technicalities of semi-analytic 
models and how they are constrained may wish to skip !|2] !|3] 
and most of !|4j and jump directly to §4.12| where we present 
two interesting predictions from our model and ^in which 
we explore the effects of varying key physical processes. 



Benson et al. (20031: Effects of thermal conduction on 



cluster cooling rates and "superwind" feedback from super- 
novae (described in further detail by[Baugh et al. 20051. 

• [Benson et al. (2004( | : Heating of galactic disks by orbit- 
ing dark matter halos. 



Nagashima et al. (2005a I: Detailed chemical enrichment 



models (incorporating delays and tracking of individual ele- 



ments; see S 2.11 1 



Bower et al. (2006[): F eedback from AGN (see g2.13[). 
Malbon et al. (2007[): Black ho le gro wth (see p.l3} 



as applied to the Baugh et al. (2005 1 — see Fanidakis et al 



(2009 1 for a similar (and more advanced) treatment of black 
hole s in the [Bower et al. (2006| ) model. 



Stringer & Benson (20071: Radially resolved structure 



of galactic disks. 

• [Font et al. (2008] ): Ram pressure stripping of cold gas 
from galactic disks (see [2.9 1. 



2.2 Executive Summary 

Having developed these treatments of various physical 
processes one-by-one, our intention is to integrate them into 
a single baseline model. In addition to the accumulation 
of many of these improvements (many of which have not 
previously been utilized simultaneously), the two major 
modifications to the Galform model introduced in this 
work are: 



2 MODEL 

In this section we provide a detailed description of our 
model. 



2.1 Starting Point 



The starting point for this discussion is [Cole et al. (2000] ) 
and we will refer to that work for details which have not 
changed in the current implementation. We choose |Cole| 
|et al. (2000[ ) as a starting point for the technical description 
of our model as it represents the last point at which the 
details of the Galform model were presented as a coherent 
whole in a single document. As noted in ^ however, the 
scientific predecessor of this work is that of [Bower et al.[ 
(20061. That paper, and several others, introduced many 
improvements relative to [Cole et al. (2000[ ), many of which 
are described in more detail here. A brief chronology of the 



development of Galform from Cole et al. (20001 to the 
present is as follows: 



Cole et al. (2000 1: Previous full description of the Gal- 



form model. 



Grasil (see [2.14.1 1. 
• [Benson et al. (2001[ ): 



Granato et al. (20001: Detailed dust modelling utilizing 



Treatment of reionization and the 



evol ution of the IGM (s ee p.lO ). 

• Bower et al. (200T| ): Treatment of heating and ejection 



of hot material from halos due to energy input (see [ 2.13 1. 
• [Benson et al. (2002b[ ): Back reaction of reionization and 



photoionizing background on galaxy formation (see [2.101 
and detailed treatment of satellite galaxy dynamics (a some- 



what different approach to this is described in [ 2.8 and [ 2.9 1 



• The removal of discrete "formation" events for dark 
matter halos (which previously occurred each time a halo 
doubled in mass and caused calculations of cooling and 
merging times to be reset). This has facilitated a major 
change in the Galform cooling model which previously 
made fundamental reference to these formation events. 

• The inclusion of arbitrarily deep levels of subhalos 
within subhalos and, as a consequence, the possibility of 
mergers between satellite galaxies. 

Aspects of the model that are essentially unchanged 
from Cole et al. (20001 are listed in [2.3 Before launch- 



ing into the detailed discussion of the model, [ ]2.4| pro- 
vides a quick overview of what has changed between [Cole[ 



et al. (2000 1 and the current implementation. In addition 
to changes to the physics of the model, the Galform code 
has been extensively optimized and made OpenMP parallel 
to permit rapid calculation of self-consistent galaxy /IGM 



evolution (see [2.101 



2.3 Unchanged Aspects 

Below we list aspects of the current implementation of Gal- 
form that are unchanged relative to that published in [Cole[ 
et al. (2000]). 



Virial Overdensities: Virial overdensities of dark mat- 



ter halos are computed as described by Cole et al. (20001, 
i.e. using the spherical top-hat collapse model for the ap- 
propriate cosmology and redshift. Given the mass and virial 
overdensity of each halo the corresponding virial radii and 
velocities are easily computed. 

• Star Formation Rate: The star formation rate in disk 
galaxies is given by 
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Mcom/t* where r* = Vdisk(Vdisk/200km s 



(1) 



where Afcoid is the mass of cold gas in the disk, Tdisk = 
J'disk/Vdisk is the dynamical time of the disk at the half mass 
radius r^isk and Vdisk is the circular velocity of the disk at 
that radius. The two parameters and control the nor- 
malization of the star formation rate and its scaling with 
galaxy circular velocity respectively. 

• Mergers/Morphological Transformation: The classifica- 
tion of merger events as minor or major follows the logic of 
Cole et al. (2000 §4.3.2). However, the rules which deter- 
mine when a burst of star formation occurs are altered to 
become: 

- Major merger? 

• Requires Msat/Mccn > /burst. 

- Minor merger? 

( A'/ccn(bulgc)/^con < ^/T^^urst 

■ Requires < and 

L A/ccn(cold) /A/ccn ^ /gas, burst • 

where Mcon and Msat are the baryonic masses of the central 
and satellite galaxies involved in the merger respectively, 



Mr. 



i(bulge) 



is the mass of the bulge component in the cen- 



tral galaxy and /burst, /gas, burst and B/Tburst are parameters 
of the model. The parameter B/Tburst is intended to inhibit 
minor merger-triggered bursts in systems that are primarily 
spheroid dominated (since we may expect that in such sys- 
tems the minor merger cannot trigger the same instabilities 
as it would in a disk dominated system and therefore be 
unable to drive inflows of gas to the central regions to fuel 
a burst). We would expect that the value of this parameter 
should be of order unity (i.e. the system should be spheroid 
dominated in order thatthe burst triggering be inhibited). 

• Spheroid Sizes: The sizes of spheroids formed through 
mergers are computed using the approach described by |Cole| 

pooo' 



et al. 



§4.4.2). 



• Calculation of Luminosities: The luminosities and mag- 
nitudes of galaxy are computed from their known stel- 
lar populations as described by Cole et al. (2000 §5.1). 



(Although note that the treatment of dust extinction has 



changed; see f 2.14.1 



2.4 Overview of Changes 

We list below the changes in the current implementation of 



Galform relative to that published in Cole et al. (20001. 
These are divided into "minor changes" , which are typically 
simple updates of fitting formulas, and "major changes", 
which are significant additions to or modifications of the 
physics and structure of the model. 



2.4-1 Minor changes 



• Dark matter halo mass function: [See [ 2.5.1| Cole et al. 
|(2000P use the [Press fc Schechter (1974] ) mass function for 
dark matter halos. In this work, we use the more recent de- 



termination of Reed et al. (20071 which is calibrated against 
N-body simulations over a wide range of masses and red- 
shifts. 



Dark matter merger trees: [See [ 2.5.2 Cole et al. (2000 1 



use a binary split algorithm utilizing halo merger rates in- 



ferred from the extended Press-Schechter formalism (Lacey 



& Cole 19931. We use an empirical modification of this al- 



gorithm proposed by Parkinson et al. (2008 1 and which pro- 
vides a much more accurate match to progenitor halo mass 
functions as measured in N-body simulations. 



Density profile of dark matter halos: [See [ 2.5.4 Cole 



et al. (20001 employed NFW (Navarro et al. 19971 density 



profiles. We instead use Einasto density profiles (Einasto 



1965 1 consistent with recent findings ( Navarro et al. 2004 



Merritt et al. 2005 Prada et al. 2006 1 



Density and 



ar momentum of halo gas: [See \ 2.6.3 



Cole et al. (2000 1 adopted a cored isothermal profile for the 



hot gas in dark matter halos and furthermore assumed a 
solid body rotation, normalizing the rotation speed to the 
total angular momentum of the gas (which was assumed to 
have the same average specific angular momentum as the 
dark matter). We choose to adopt the density and angular 
momentum distributions measured in hydrodynamical sim- 
ulations by Sharma k, Steinmetz (2005| . 



Dynamical friction timescales: [See [ 2.8.5 Cole et al 



(20001 estimated dynamical friction timescales using the 



expression derived by Lacey & Cole (19931 for isothermal 



dark matter halos and the distribution of orbital parame- 



ters found by Tormen (19971. In this work, we adopt the 



fitting formula of Jiang et al. (2008 1 to compute dynamical 
friction timescales and the orbital parameter distribution of 



Benson (2005 1 



Disk stability: We retain the same test of disk stability 



as did Cole et al. (2000 1 and similarly assume that unstable 



disks undergo bursts of star formation resulting in the forma- 
tion of a spheroicj^ One slight difference is that we assume 
that the instability occurs at the largest radius for which the 
disk is deemed to be unstable rather than at the rotational 
support radius as [Cole et al. (2000[ ) assumed. This prevents 
galaxies with very low angular momenta from contracting 
to extremely small sizes (and thereby becoming very highly 
self-gravitating and unstable) before the stability criterion is 
tested. Additionally, we allow for different stability thresh- 
olds for gaseous and stellar disks. We employ the stability 
criterion of [Efstathiou et al. (1982[ ) such that disks require 

Vd 

(GMdAR^ ^ 

to be stable, where Vd is the disk rotation speed at the half- 
mass radius, Md is the disk mass and Rs is the disk radial 



scale length. Efstathiou et al. (1982 \ found a value of ed 



1.1 was applicable for purely stellar disks. |Christodoulou| 



et al. (1995 1 demonstrate that an equivalent result for 
gaseous disks gives td.gas = 0.9. We choose to make td.gas a 
free parameter of the model and enforce td,* ~ ed,gas + 0.2. 
For disks containing a mixture of stars and gas we linearly in- 
terpolate between td,* and td.gas using the gas fraction as the 
interpolating variable. As has been recently pointed out by 
[Athanassoula (20081, this treatement of the process of disk 
destabilization, similar to that in other semi-analytic mod- 



els, is dramatically oversimplified. As Athanassoula (2008 \ 



^ While the implementation of this physical process is unchanged, 
[Cole et al. (2000] l actually ignored this process in their fiducial 
model, while we include it in our work. 
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also describes, a more realistic model would need both a 
much more careful assessment of the disk stability and a 
consideration of the process of bar formation. This currently 
remains beyond the ability of our model to address, although 
it should clearly be a priority area in which semi-analytic 
models should strive to improve. In Galform we can con- 
sider an alternative disk instability treatment in which dur- 
ing an instability event only just enough mass is transferred 
from the disk to the spheroid component to restabilize the 
disk. While this does not explore the full range of uncer- 
tainties arising from the treatment of this process, it gives 
at least some idea of how significant they may be. We find 
that the net result of switching to the alternative treatment 
of instabilities is to slightly increase the number of bulgeless 
galaxies at all luminosities, with a corresponding decrease 
in the numbers of intermediate and pure spheroid galaxies. 
The changes, however, do not alter the qualitative trends 
of morphological mix with luminosity nor global properties 
of galaxies such as sizes and luminosity functions at z = 0. 
At higher redshifts (e.g. 2 > 5), the change is more signifi- 
cant, with a reduction in star formation rate by a factor of 
2-3 resulting from the lowered frequency of bursts of star 
formation. This change could be offset by adjustments in 
other parameters, but demonstrates the need for a refined 
understanding and modelling of the disk instability process 
in semi-analytic models. 

• Sizes of galaxies: [See [ 2.7 . Sizes of disks and spheroids 
are determined as described by [Cole et al. (2000| ), although 
the equilibrium is solved for in the potential corresponding 
to an Einasto density profiles (used throughout this work) 
rather than the NFW profiles assumed by [Cole et al. (2000] ) 
and adiabatic contraction is computed using the methods of 



Gnedin et al. (20041 rather than that of Blumenthal et al, 
(1986| ). 



While we class the above as minor changes, the effects of 
some of these changes can be significant in the sense that re- 
verting to the previous implementation would change some 
model predictions by an amount comparable to or greater 
than the uncertainties in the relevant observational data. 
However, none of these modifications lead to fundamental 
changes in the behaviour of the model and their effects could 
all be counteracted by small adjustments in model param- 
eters. This is why we classify them as "minor" and do not 
explore their consequences in any greater detail. 



2.4-2 Major changes 



Spins of dark matter halo: [See [ 2.5.5 In Cole et al 



(20001 spins of dark matter halos were assigned randomly 
by drawing from the distribution of |Cole fc Lacey (1996[ ). In 
this work, we implement an updated version of the approach 



described by Vitvitska et al. (2002 I to produce spins corre- 
lated with the merging history of the halo and consistent 



(20071 



with the distribution measured by Bett et al. 

• Removal of discrete formation events: [See [2.5.3 The 
discrete "formation" events (associated with mass dou- 



blings) in merger trees which Cole et al. (2000) ) utilized to 
reset cooling and merging calculations are no more. Instead, 
cooling, merging and other processes related to the merger 
tree evolve smoothly as the tree grows. 

• Cooling model: [See [ 2.6 The cooling model has been 



revised to remove the dependence on halo formation events, 
allow for gradual recooling of gas ejected by feedback and 
accounts for cooling due to molecular hydrogen and Comp- 
ton cooling and for heating from a photon background. 
• Ram pressure and tidal stripping [See [2.9 Ram pres- 



sure and tidal stripping of both hot halo gas and stars and 
interstellar medium (ISM) gas in galaxies are now accounted 
for. 



• IGM interaction [See [2.10 Galaxy formation is solved 
simultaneously with the evolution of the intergalactic 
medium in a self-consistent way: emission from galaxies and 
AGN ionize and heat the IGM which in turn suppresses the 
formation of future generations of galaxies. 

• Full hierarchy of subhalos [See [2.8 All levels of the 
substructure hierarchy (i.e. subhalos, sub-subhalos, sub-sub- 
subhalos. . . ) are included in calculations of merging. This 
allows for satellite-satellite mergers. 

• Non-instantaneous recycling [See [ 2.11| The instanta- 
neous recycling approximation for mass loss, chemical en- 
richment and feedback has been dropped and the full time 
and metallicity-dependences included. All models presented 
in this work utilize fully non-instantaneous recycling, metal 
production and supernovae feedback. 



2.5 Dark Matter Halos 

2.5.1 Mass Function 

We assume that the masses of dark matter halos at any given 
redshift are distributed according to the mass function found 
by |Reed et al. (2007| |. Specifically, the mass function is given 
by 



An 



dluMv 



2 SloPcrit 
TT 



diner 



x[l + 1.047(a; 



X exp 



1 2 
'2" 



dluM 
^'') + O.6G1 + 0.4G2]^'a; 

2p 



0.0325 



UJ 



{n,s + 3)2 



(3) 



where dn/d In is the number of halos with virial mass 
per unit volume per unit logarithmic interval in Mv, a{M) is 
the fractional mass root-variance in the linear density field 
in top-hat spheres containing, on average, mass M, 5c{z) 
is the critical overdensity for spherical top-hat collapse at 



redshift z ( Eke et al. 1996 \ 



ncB 



dlna 
dluM 

Sciz) 



Gi 



exp 



(4) 
(5) 

(6) 

(7) 
(8) 

A' = .310, ca = 0.764 a nd p = 0.3 as in eqns. (11) and 
(12) of |Reed et al. (2007| p| The mass variance, cr^(M), is 



G2 = exp 




With minor corrections to the published version (Reed, private 
communication) . 
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computed using the cold dark matter transfer function of 



Eisenstein & Hu (1999 1 together with a scale- free primordial 
power spectrum of slope ris and normalization erg. 

When constructing samples of dark matter halos we 
compute the number of halos, Ai'haio, expected in some vol- 
ume V of the Universe within a logarithmic mass interval, 
A In My, according to this mass function, requiring that the 
number of halos in the interval never exceeds A'^max and is 
never less than A^'min to ensure a fair sample. We then gener- 
ate halo masses at random using a Sobol' sequence ( [Sobol'l 
1967 1 drawn from a distribution which produces, on average. 



A'haio halos in each interval. This ensures a quasi-random, 
fair sampling of halos of all masses with no quantization of 
halo mass and with sub-Poissonian fluctuations in the num- 
ber of halos in any mass interval. 



2.5.2 Merger Trees 

Dark matter halo merger trees, which describe the hierar- 
chical growth of structure in a cold dark matter universe, 
form the backbone of our model within which the process 
of galaxy formation proceeds. Merger trees are either con- 
structed through a variant of the extended Press-Schechter 
Monte Carlo methodology, or are extracted from N-body 
simulations. 

When constructing trees using Monte Carlo methods, 
we employ the merger tree algorithm described by |Parkm^ 
[son et al. (2008"1 ) which is itself an empirical modification 
of that described by Cole et al. (20001. We adopt the pa- 
rameters (Go, 7i, 72) = (0.57,0.38,-0.01] that Parkinson 



et al. (2008 1 found provided the best fil|^ to the statistics 
of halo progenitor masses measured from the Millennium 
Simulation by Cole et al. (2008)) . We typically use a mass 
resolution (i.e. the lowest mass halo which we trace in our 
trees) of 5 x W^h~^ Mq, which is sufficient to achieve re- 
solved galaxy properties for all of the calculations consid- 
ered in this work. An exception is when we consider Local 



Group satellites (see S 4.10 1, for which we instead use a mass 
resolution of IO^/i'^Mq. Figure [l] shows the resulting dark 
matter halo mass functions at several different redshifts and 
demonstrates that they are in good agreement with that 
expected from eqn. ([3|. 




Figure 1. The dark matter halo mass function is shown at a 
number of different redshifts. Solid lines indicate the mass func- 
tion expected from eqn. |3| while points with error bars indicate 
the mass function constructed using merger trees from our model. 
The trees in question were initiated at z = and grown back to 
higher redshifts using the methods of |Parkinson et al. 



formation event). Additionally, any gas ejected by feedback 
was only allowed to begin recooling after a formation event, 
and any satellite halos that had not yet merged with the 
central galaxy of their host halo were assumed to have their 
orbits randomized by the formation event and consequently 
their merger timescales were reset. 

While computationally useful, these formation events 
lack any solid physical basis. As such, we have excised them 
from our current implementation of Galform. Halo prop- 
erties (virial velocity and spin) now change at each timestep 
in response to mass accretion. Additionally, the cooling and 
merging calculations no longer make use of formation events 
(see j]2.6|and j^2.8| respectively) . 



2.5.3 (Lack of) Halo Formation Events 



Cole et al. (2000 1 identified certain halos in each dark mat- 
ter merger tree as being newly formed. "Formation" in this 
case corresponded to the point where a halo had doubled in 
mass since the previous formation event. The characteristic 
circular velocity and spin of halos was held fixed in between 
formation events and the time available for hot gas in a halo 
to cool was measured from the most recent formation event 
(such that the cooling radius was reduced to zero at each 



^ [Benson (2008l l found an alternative set of parameters which 
provided a better match to the evolution of the overall halo mass 
function but performed slightly less well (although still quite well) 
for the progenitor halo mass functions. We have chosen to use 
the parameters of [Parkinson et al. (2008l l as for the properties 
of galaxies we wish to get the progenitor masses as correct as 
possible. 



2.5.4 Density Profiles 



200. 



Recent N-body studies (Navarro et al. 200^ VIerritt et al 
l^rada et al. 200^ 1 indicate that the density profiles of 



dark matter halos in CDM cosmologies are better described 



profile, 



by the Einasto profile ( Eineisto 1965 1 than the NFW profile 



(Navarro et al. 19971. As such, we use the Einasto density 



p(r) 



P-2 exp 



r 

r-2 



1 



(9) 



where r_2 is a characteristic radius at which the logarithmic 
slope of the density profile equals —2 and a is a parameter 
which controls how rapidly the logarithmic slope varies with 
radius. To fix the value of a we adopt the fitting formula of 



Gao et al. (20081, truncated so that a never exceeds 0.3, 



0.155 + 0.0095;/^ 
0.3 



ifi^ < 3.907 
■du> 3.907, 



(10) 
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where u = Sc{a) /a{M) which is a good match to halos in the 
Millennium SimulatiorQ The value of r_2 for each halo is 
determined from the known virial radius, r^, and the concen- 
tration, c_2 = rv/r_2. Concentrations are computed using 



the method of Navarro et al. (19971 but with the best-fit 



parameters found by Gao et al. (20081. 

Various integrals over the density and mass distribution 
are needed to compute the enclosed mass, angular momen- 
tum, velocity dispersion, gravitational energy and so on of 
the Einasto profile. Some of these may be expressed analyti- 



cally in terms of incomplete gamma functions ( Cardone et al 



2005). Expressions for the mass and gravitational potential 
are provided by [Cardone et al. ( 2005 ). One other integral, 
the angular momentum of material interior to some radius, 
can also be found analytically: 



J(r) 



Jo 

2yr 4 



'(3 + <»rot)„/^'NJ / 



p(r )dr 



(f) 



xr 



4 + Q,ot 2(r/r-_2)" 



(11) 



where the specific angular momentum at radius r is assumed 
to be rVrotij- / r^)"'^°^ and T is the lower incomplete gamma 
function. Other integrals (e.g. gravitational energy) are com- 
puted numerically as needed. 



2.5.5 Angular momentum 



As first suggested by Hoyle (19491, and developed further 
by [Doroshkevich (1970| ), [Peebles (1969t and [White (1984] ), 
the angular momenta of dark matter halos arises from tidal 
torques from surrounding large scale structure and is usually 
characterized by the dimensionless spin parameter. 



A ; 



Jvl-Ev 



1/2 



CM, 



5/2 



(12) 



where Jv is the angular momentum of the halo and its 
energy (gravitational plus kinetic). The distribution of A 
has been measured numerous times from N-body simula- 



tions (Barnes & Efstathiou 1987 Efstathiou et al. 198J iVar 
ren et al. 199S Ziole fc Lacey 199(^|Lemson &: KaufFmann 1999 1 



and found to be reasonably well approximated by a log- 
normal distribution. More recent estimates by [Belt et al.| 
(20071 using the Millennium Simulation show a somewhat 



different form for this distribution: 

Ao. 



exp 



(13) 



where a\ = 2.509 and Ao = 0.04326 are parameters. 

[Cole et al. (2000[ ) assigned spins to dark matter halos 
by drawing them at random from the distribution of |Cole| 
& Lacey (19961. This approach has the disadvantage that 
spin is not influenced by the merging history of a given dark 
matter halo and, furthermore, spin can vary dramatically 



6 [Gao et al. (2008[ l were not able to probe the behaviour of a in 
tlie very liigh v regime. Extrapolating their formula to f > 4 is 
not justified and we instead choose to truncate it at a maximum 
of a = 0.3. 



from one timestep to the next even if a halo experiences 
no (or only very minor) merging. This was not a problem 
for Cole et al. (20001, who made use of the spin of each 



newly formed halo, ignoring any variation between forma- 
tion event^ However, in our case, such behaviour would 
be problematic. We therefore revisit an idea first suggested 



by [Vitvitska et al.[ ([2002| see also [Mailer et al. 2002]. They 



followed the contribution to the angular momentum of each 
halo from its progenitor halos (which carry angular momen- 
tum in both their internal spin and orbit). Note that the 
angular momentum still arises via tidal torques (which are 
responsible for the orbital angular momenta of merging ha- 
los). 

Halos in the merger tree which have no progenitors are 
assigned a spin by drawing at random from the distribution 
of Bett et al. (20071. For halos with progenitors, we proceed 
cis follows: 

(i) Compute the internal angular momenta of all progen- 



itor halos using their previously assigned spin and eqn. ( 12 1 



(ii) Select orbital parameters (specifically the orbital an- 
gular momentum) for each merging pair of progenitors by 
drawing at random from the distribution found by [Benson] 
1120051); 

(iii) Sum the internal and orbital angular momenta of all 
progenitors assuming no correlations between the directions 
of these vector^ 

(iv) Determine the spin parameter of the new halo from 
this summed angular momentum and eqn. (12 I. 



Benson (2005 1 report orbital velocities for merging halos and 



give expressions for the angular momenta of those orbits as- 
suming point mass halos. While this will be a reasonable 
approximation for high mass ratio mergers it will fail for 
mergers of comparable mass halos. In addition, halo merg- 
ers may not necessarily conserve angular momentum in the 
sense that some material, plausibly with the highest specific 
angular momentum, may be thrown out during the merging 
event leaving the final halo with a lower angular momentum. 
To empirically account for these two factors we divide the 
orbital angular momentum by a factor of /2 = 1 -|- M2/M1 
(where M2 < Mi are the masses of the dark matter halos). 
We find that this empirical factor leads to good agreement 
with the measured N-body spin distribution, but could be 
justified more rigorously by measuring the angular momen- 
tum (accounting for finite size effects) of the progenitor and 
remnant halos in N-body mergers. 

To test the validity of this approach we generated 51625 
Monte Carlo realizations of merger trees drawn from a halo 
mass function consistent with that of the Millennium Sim- 
ulation and with a range of masses consistent with that for 



which Bett et al. (2007) were able to measure spin param- 
eters and applied the above procedure. Figure [2| shows the 
results of this test. We find remarkably good agreement be- 



tween the distribution of spin measured by Bett et al. (2007 ) 



and the results of our Monte Carlo model. It should be noted 

^ As it seems reasonable to assume that the spins of a halo at 
two successive formation events, i.e. separated by a factor of two 
in halo mass, would be only weakly correlated. 

Additionally, we are assuming that mass accretion below the 
resolution of the merger tree contributes the same mean specific 
angular momentum as accretion above the resolution. 
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B07 TREE halos 
B07 SO halos 
B07 For halos 
This work 




Figure 2. The distribution of dark matter halo spin parame- 
ters. Black lines show measurements of this distribution from the 
Millennium N-body simulation |Bett et aril|2007| B07 ), for three 
different group finding algorithms. |Bett et al. (2007| note that 
the "tree" halos give the most accurate determination of the 
spin distribution. The red line shows the results of the Monte 
Carlo model described in this work, using 51625 Monte Carlo re- 
alizations of merger trees spanning a range of masses identical to 
that used by |Bett et al. (2007| |. 



that our assumption of no correlation between the various 
angular momenta vectors of progenitor halos is not correct. 
However, [Benson (2005| ) shows that any such correlations 
are weak. Therefore, given the success of a model with no 
correlations, we choose to ignore them. 

Our results are in good agreement with previous at- 
tempts to model the halo spin distribution in this way. 



Mailer et al. (2002 1 found good agreement with N-body re- 



sults using the same principles, although they found that 
introducing some correlation between the directions of spin 
and orbital angular momenta improved their fit. |Vitvitsk"a| 
et al. (2002 1 also found generally good agreement with N- 



body simulations using orbital parameters of halos drawn 
from an N-body simulation. Both of these earlier calcula- 
tions relied on much less well calibrated orbital parameter 
distributions for merging halos and the simulations to which 
they compared their results had significantly poorer statis- 
tics than the Millennium simulation. Our results confirm 
that this approach to calculating halo spins from a merger 
history still works extremely well even when confronted with 
the latest high-precision measures of the spin distribution. 



2.6 Cooling Model 



The cooling model described by Cole et al. (2000 \ deter- 



mines the mass of gas able to cool in any timestep by fol- 
lowing the propagation of the cooling radius in a notional 
hot gas density profile]^ which is fixed when a halo is fiagged 



as "forming" and is only updated when the halo undergoes 
another formation event. The mass of gas able to cool in any 
given timestep is equal to the mass of gas in this notional 
profile between the cooling radius at the present step and 
that at the previous step. The cooling time is assumed to 
be the time since the formation event of the halo. Any gas 
which is reheated into or accreted by the halo is ignored un- 
til the next formation event, at which point it is added to 
the hot gas profile of the newly formed halo. The notional 
profile is constructed using the properties (e.g. scale radius, 
virial temperature etc.) of the halo at the formation event 
and retains a fixed metallicity throughout, corresponding to 
the metallicity of the hot gas in the halo at the formation 
event. 

In this work we implement a new cooling model. We do 
away with the arbitrary "formation" events and instead use 
a continuously updating estimate of cooling time and halo 
properties. For the purposes of this calculation we define the 
following quantities: 

• Afhot: The current mass of hot (i.e. as yet uncooled) gas 
remaining in the notional profile; 

• Afcooicd: The mass of gas which has cooled out of the 
notional profile into the galaxy phase; 

• Mrehcatod: The mass of gas which has been reheated (by 
supernovae feedback) but has yet to be reincorporated back 
into the hot gas component. 

• Mcjoctcd: The mass of gas which has been ejected beyond 
the virial radius of this halo, but which may later reaccrete 
into other, more massive halos. 

The notional profile always contains a mass Mtotai ~ 
Mhot + Afcooied + A/rehoated- The properties (density normal- 



ization, core radius) are reset, as described in [ 2.6.3 at each 



timestep. The previous infall radius (i.e. the radius within 
which gas was allowed to infall and accrete onto the galaxy) 
is computed by finding the radius which encloses a mass 
Afcooicd -I- Mrchoatcd (i-c. the mass previously removed from 
the hot component) in the current notional profile. 

We aim to compute a time available for cooling for the 
halo, tavaii, from which we can compute a cooling radius 
in the usual way (i.e. by finding the radius in the notional 
profile at which tcooi = iavaii)- In Cole et al. (20001 the time 
available for cooling is simply set to the time since the last 
formation event of the halo. 

At any time, the rate of cooling per particle is just 
A(r, Z, riH, F^)nH where A(r, Z, nn, Fv) is the cooling func- 
tion, and nu the number density of hydrogen, Z a vector of 
metallicity (such that the i*'^ component of Z is the abun- 
dance by mass of the i^^ element) and the spectrum of 
background radiation. The total cooling luminosity is then 
found by multiplying by the number of particles, A^, in some 
volume V that we want to consider. If we take this volume to 
be the entire halo then N = Mtota.1 / fJ-mn- If we integrate this 
luminosity over time, we find the total energy lost through 
cooling. The total thermal energy in our volume V is just 
3NkBT/2. The gas will have completely cooled once the en- 
ergy lost via cooling equals the original thermal energy, i.e.: 



^ We refer to this as a "notional" profile since it is taken to 
represent the profile before any cooling can occur. Once some 



cooling occurs presumably the actual profile adjusts in some way 
to respond to this and so will no longer look like the notional 
profile, oven outside of the cooling radius. 
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3NkBn/2 = / A{t')nnNdt', 



where for brevity we write A(t) 
A[T^{t),Z{t),nH{t),F^{t)]. We can write this as 



^cool ^avail ' 



where 



2A(t)nH 
is the usual cooling time and 



(14) 



(15) 



(16) 



^ avail 



J^A{t')nu{t')Ndt' 
Ait)nn{t)N 



(17) 



is the time available for cooling. We can re-write this as 
_ /J[Tv(t')A^Acooi(t')]dt' 



[rv(t)Ar/tcooi(t)] 



(18) 



In the case of a static halo, where Tv, Z, and A'' are 
independent of time, tavaii reduces to the time since the halo 
came into existence as we might expect. For a non-static halo 
the above makes more physical sense. For example, consider 
a halo which is below the 10*K cooling threshold from time 
t = to time t — t^, and then moves above that threshold 
(with fixed properties after this time). Since fcooi ~ oo (i.e. 
A{t) = 0) before in this case we find that favaii = t — 
as expected. Note that since the number of particles, A'^, 



appears in both the numerator and denominator of eqn. ( 18 1 
we can, in practice, replace A' by Mtotai without changing 
the resulting time. 

The cooling time in the above must be computed for 
a specific value of the density. We choose to use the cool- 
ing time at the mean density of the notional profile at each 
timestep. This implicitly assumes that the density of each 
mass element of gas in the notional profile has the same 
time dependence as the mean density of the profile, i.e. that 
the profile evolves in a self-similar way and that A{t) is in- 
dependent of nu (which will only be true in the collisional 
ionization limit). This may not be true in general, but serves 
as an approximation allowing us to describe the cooling of 
the entire halo with just a single integrap°1 

Having computed the time available for cooling we 
solve for the cooling radius in the notional profile at which 
tcooi(fcooi) = iavaii (as described in { 2.6.41. We also estimate 
the largest radius from whch gas has had sufficient time to 
freefall to the halo centre (as described in S 2.6.5 1. The cur- 



rent infall radius is taken to be the smaller of the cooling 
and freefall radii. Any mass between the current infall ra- 
dius and that at the previous timestep is allowed to infall 
onto the galaxy during the current timestep — that is, it is 
transferred from Mhot to Mcooicd- 

One refinement which must be introduced is to limit the 
integral 



fee / [rv(t')A-Acooi(t')]di', 



(19) 



^" A more elaborate model could compute a separate integral 
for each shell of gas, following the evolution of its density as a 
function of time as the profile evolves due to continued infall and 
cooling. 



SO that the total radiated energy cannot exceed the total 
thermal energy of the halo. This limit is given by 

fmax = Ik^T.im P'°'f , , (20) 
2 Ptotal('"v) 

where ptotai is the mean density of the notional profile and 
ptotaiij'v) is the density of the notional profile at the virial 
radius. For the entire halo (out to the virial radius) to cool 
takes longer than for gas at the mean density of the halo to 
cool, by a factor of Ptotai/ptotai(''v)- This is the origin of the 
ratio of densities in eqn. ( 20 1 . 

We must then consider two additional effects: accretion 



(S 2.6.1 1 and reheating (i 2.6.2 1. The cooling model is then 



fully specified once we specify the distribution of gas in the 



notional profile (S 2.6.3 1, determine a cooling radius (S 2.6.4 1 



and freefall radius (i 2.6.5 1, and consider how to compute 
the angular momentum of the infalling gas (f 2.6.6 1. 



2. 6. 1 Accretion 

When a halo accretes another halo, we merge their notional 
gas profiles. Since the integral, £ = J {NT^ /tcooi)dt, that we 
are computing is the total energy lost we simply add £ from 
the accreted halo to that of the halo it accretes into. This 
gives the total energy lost from the combined notional pro- 
file. However, we must consider the fact that only a fraction 
Mhot /Mtotai of the gas from the accreted halo is added to 
the hot gas reservoir of the combined halo (the mass Mcooied 
from the accreted halo becomes the satellite galaxy while the 
mass Mrcheatod is added to the reheated reservoir of the new 
halo to await reincorporation into the hot component; see 
S 2.6.2 I. We simply multiply the integral £ of the accreted 



halo by this fraction before adding it to the new halo. 

Figure [3] compares the mean cooled gas fractions in ha- 
los of different masses computed using the cooling model 
described here (green lines) and two previous cooling mod- 
els used in Galform: that of Cole et al. (2000 red lines) 



and that of Bower et al. (2006 blue lines). The only sig- 



nificant difference between the cooling implementations of 
Cole et al. (2000|) and [Bower et al. (2006) is that [Bower 



et al. (2006 1 allow reheated gas to gradually return to the 



hot component (and so be available for re-cooling) at each 
timestep (in the same manner as in the present work), while 
Cole et al. (2000 1 simply accumulated this reheated gas and 
returned it all to the hot component only at the next halo 
formation event (i.e. after a halo mass doubling). No star or 
black hole formation was included in these calculations, so 
consequently there is no reheating of gas, expulsion of gas 
from the halo or metal enrichment. Additionally, no galaxy 
merging was allowed. The thick lines show the total cooled 
fraction in all branches of the merger trees, while the thin 
lines show the cooled fraction in the main branch of the 
tree£3 



The cooling model utilized by Bower et al. (2006 \ was 



We define the main branch of the merger tree as the set of pro- 
genitor halos found by starting from the final halo and repeatedly 
stepping back to the most massive progenitor of the current halo 
at each time step. It should be noted that definition is not unique, 
and can depend on the time resolution of the merger tree. It can 
also result in situations where the main branch does not corre- 
spond to the most massive progenitor halo at a given timestep. 
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similar to that of Cole et al. (2000 1 except that it allowed 



accreted and reheated gas to rejoin the hot gas reservoir in 
a continuous manner rather than only at each halo forma- 
tion event. Additionally, it used the current properties of 
the halo (e.g. virial temperature) to compute cooling rates 
rather than the properties of the halo at the previous for- 



mation event. As such, the Bower et al. (20061 model con- 



tains many features of the current cooling model, but retains 
the fundamental division of the merger tree into discrete 
branches as in the Cole et al. (2000 1 model. 

We find that, in general, the cooling model described 
here predicts a total cooled fraction very close to that pre- 



dicted by the cooling model of Bower et al. (20061, the ex- 
ception being at very early times in low mass halos where it 
gives a slightly lower value. The difference of course is that 
the new model does not contain artificial resets in the cool- 
ing calculation which, although they make little difference 
to this statistic, have a strong influence on, for example, 
calculations of the angular momentum of cooling gas. Both 
of these models predict somewhat more total cooled mass 



than the Cole et al. (20001 model. This is due entirely to 



the allowance of accreted gas to begin cooling immediately. 

If we consider the cooled fraction in the main branch 
of each tree (i.e. the mass in what will become the central 
galaxy in the final halo) we see rather different behaviour. 
At early times, the new model tracks the Bower et al. (2006 1 
model. At late times, however, the |Bower et al. (2006[ ) model 
shows a much lower cooling rate while the new model tracks 



the cooled fraction in the Cole et al. (2000 1 model quite 



closely. This occurs in massive halos where, in the |Bower| 
et al. (2006 1 model the use of the current halo properties 
to determine cooling rates results in ever increasing cooling 
times as the virial temperature of the halo increases and the 



halo density (and hence hot gas density) decline. The Cole 



et al. (2000 1 model is less susceptible to this as it computes 



halo properties based on the halo at formation. The new 
cooling model produces results comparable to the |Cole et al.| 
(20001 model since, while it utilizes the present properties 
of the halo just as does the [Bower et al. (2006[ ) model, it 
retains a memory of the early properties of the halo. 



2.6.2 Reheating 



When gas is reheated (via feedback; S 2.12 1 we assume that 



it is heated to the virial temperature of the current halo 
(i.e. the host halo for satellite galaxies) and is placed into a 
reservoir Mreheated. Mass is moved from this reservoir back 
into the hot gas reservoir on a timescale of order the halo 
dynamical time, Tdyn. Specifically, mass is returned to the 
hot phase at a rate 



Mr, 



Tdyn 



(21) 



during each timestep. This effectively undoes the cooling 
energy loss which caused this gas to cool previously. The 
energy integral £ is therefore modified by subtracting from 
it an amount AA'^rchcatodTv where AA'^rchoatcd is the number 
of particles reheated. 

Similarly, the notional profile is allowed to "forget" 
about any cooled gas on a timescale of order the dynam- 
ical time (i.e. we assume that the notional profile adjusts to 




Figure 3. The mean cooled gas fractions in the merger trees 
of halos with masses lO", lO^^, lo", IQl" and IO^^H-^Mq at 
2 = are shown by coloured lines. Green lines show results from 
the cooling model described in this work while red lines indicate 
th e model of|Cole et al. (2000| and blue lines the cooling model 
of [Bower et al. (2006[ |. Solid lines show the total cooled fraction 
in all branches of the merger trees, while dashed lines show the 
cooled fraction in the main branch of the trees. For the purposes 
of this figure, no star or black hole formation was included in 
these calculations, so consequently there is no reheating of gas, 
expulsion of gas from the halo or metal enrichment. Additionally, 
no galaxy merging was allowed. As such, the differences between 
models arises purely from their different implementations of cool- 
ing. 



the loss of this gas). This is implemented by removing mass 
from the cooled reservoir at a rate 

-^'-^coolcd 



•t remove 



(22) 



2.6.3 Hot Gas Distribution 

The hot gas is assumed to be distributed in a notional profile 
with a run of density consistent with that found in hydro- 



dynamical simulations ( Sharma & Steinmetz 200E itringer & 
[Benson 2007[ ). [Sharma fc Steinmetz (2005[ ) performed non- 
radiative cosmological spectral energy distribution (SPH) 
simulations and studied the properties of the hot gas in dark 
matter halos. These simulations are therefore well suited to 
our purposes since they relate to the notional profile which 
is defined to be that in the absence of any cooling. The gas 
density profiles found by [Sharma &: Steinmetz (2005[ ) are 
well described by the expression: 
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p(r) 



(r- 



(23) 



where rcore is a characteristic core radius for the profile. We 
choose to set rcore = icorerv whcre Ocore is a parameter whose 
value is the same for all halos at all redshifts. The simula- 
tions suggest that Ocorc ~ 0.05 (Stringer & Benson 20071, 
but we will treat ttcorc as a free parameter to be constrained 
by observational data. The density profile is normalized such 
that 



p(r)47rr dr = Mtotai, 



(24) 



and the hot gas is assumed to be isothermal at the virial 
temperature 



_ 1 fjmH 2 
^"'2 k 



(25) 



with a metallicity equal to Z = M^^hot/Afhot. Initially, 
Afz.hot = but can become non-zero due to metal produc- 
tion and outflows as a result of star formation and feedback. 



2,.6.^ Cooling Radius 



Given the time available for cooling from eqn. ( 18 1 the cool- 
ing radius is found by solving 



f (ntot/".H)fcBrv 
nH(rcooi)A(t) 



(26) 



where ntot is the total number density of the atoms in the 
gas. Due to the dependence of A(t) on density when a pho- 



toionizing background is present (see S5.ll this equation 
must be solved numerically. 



2.6.5 Freefall Radius 

To compute the mass of gas which can actually reach the 
centre of a halo potential well at any given time we require 
that not only has the gas had time to cool but also that it 
has had time to freefall to the centre of the halo starting 
from zero velocity at its initial radius. To estimate the max- 
imum radius from which cold gas could have reached the 
halo centre through freefall we proceed as follows. We com- 
pute a time available for freefall in the halo, tavaii,ff, using 



eqn. (18 1, but limit the integral £ (defined in eqn. 191 such 



that the time available can not exceed the freefall time at 
the virial radius. We then solve the freefall equation 

tavail,ff, (27) 



^2[$(r') - -l'(rff)] 

where <E>(r-) is the gravitational potential of the halo, for the 
radius rg at which the freefall time equals the time available. 
Only gas within the minimum of the cooling and freefall radii 
at each timestep is allowed to reach the centre of the halo 
and become part of the forming galaxy. 



Jhot: The total angular momentum in the Mhot reservoir 
of the notional profile; 

Jcooied: The total angular momentum in the Afcooied reser- 
voir of the notional profile; 

J'rchcatcd: The total angular momentum in the Mrohcated 
reservoir of the notional profile; 

inew: The specific angular momentum which newly ac- 
creted material must have in order to produce the correct 
change in angular momentum for this halcf^ 

(^cooled and Jyc heated are initialized to zero at the start of the 
calculation. Jhot is initialized by assuming that any material 
accreted below the resolution of the merger tree arrives with 
the mean specific angular momentum of the halo. Angular 
momentum is then tracked using the following method: 

(i) At the start of a time step, all three angular momen- 
tum reservoirs from the most massive progenitor halo are 
added to those of the current halo. 

(ii) We assume that the specific angular momentum of the 
gas halo is distributed according to the results of |Sharma fc| 
[Stcinmetz (2005 1 such that the differential distribution of 
specific angular momentum, j, is given by 



1 dM 



-3 Ha 



(28) 



where V is the gamma function, M is the total mass of gas, 
jd = jtot I o. and jtot is the mean specific angular momentum 
of the gas. The parameter aj is chosen to be 0.89, consistent 
with the median value found by'Sharma fc Steinmetz (2005| ) 
in simulated halos. The fraction of mass with specific angular 
momentum less than j is then given by 



7 I , — 

Jd 



(29) 



where 7 is the incomplete gamma function. Once the mass 
of gas cooling in any given timestep is known the above 
allows the angular momentum of that gas to be found. This 
amount is added to the Jcooied reservoir. 

(iii) If Jroheated > thcu au augular momentum 

» T J JrehcatedCKreheat /'^dyn if Q^rcheat ^ '^dyn /on\ 

AJhot = < ^ -f A/ > 

I ^/reheated U Q^rcheat ^ '^dyn 

is transferred to back to the hot phase, consistent with the 
fraction of mass returned to the hot phase (see |2.6.2| . 

(iv) When a halo becomes a satellite of a larger halo, Jhot 
of the larger halo is increased by an amount, jnewAJhot.sat. 
This accounts for the orbital angular momentum of the gas 
in the satellite halo assuming that, on average, satellites 
have specific angular momentum of jnew. We do the same 
for Jreheated, assumiug that the Mrehoatod reservoir of the 
satellite arrives with the same specific angular momentum. 

(v) When gas is ejected from a galaxy disk to join the 
reheated reservoir it is ejected with the mean specific angular 
momentum of the disk. Gas ejected during a starburst is also 



2,. 6. 6 Angular Momentum 

The angular momentum of gas in the notional halo is tracked 
using a similar approach as for the mass. We define the fol- 
lowing quantities: 



The angular momentum of a halo differs from that of its main 
progenitor due to an increase in mass, change in virial radius and 
change in spin parameter, jnew is computed by finding the differ- 
ence in the angular momentum of a halo and its main progenitor 
and dividing by their mass difference. Note that this quantity can 
therefore be negative. 
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assumed to be ejected with the mean pseudo-specific angular 
momentunp^ of the bulge. 

Because jncw can be negative on occasion it is possible 
that Jhot < can occur. This, in turn, can lead to a galaxy 
disk with a negative angular momentum. We do not consider 
this to be a fundamental problem due to the vector nature 
of angular momentum. When computing disk sizes we sim- 
ply consider the magnitude of the disk angular momentum, 
ignoring the sign. 



2. 6. 7 Cooling/Heating Rates of Hot Gas in Halos 

The cooling model described above requires knowledge of 
the cooling function, A(r, Z, nn ,Fu). Given a gas metallicity 
and density and the spectrum of the ionizing background we 
can compute cooling and heating rates for gas in dark mat- 
ter halos. Calculations were performed with version 08.00 of 
Cloudy, last described by ,Ferland et al. (IQQSJ . In practice, 
we compute cooling/heating rates as a function of tempera- 
ture, density and metallicity using the self-consistently com- 



puted photon background (S 2.10 1 after each timestep. The 
rates are computed on a grid which is then interpolated on 
to find the cooling/heating rate for any given halo. 

Chemical abundances are assumed to behave as follows: 

• Z = : "zero" metallicity corresponding to the "pri- 
mordial" abundance ratios as used by Cloudy version 08.00 
(see the Hazy documentation of Cloudy for details). 

• [Fe/H]< —1 : "primordial" abundance ratios from 



Sutherland fc Dopita (1993| ); 

• [Fe/H]> 1 : Solar abundance ratios as used by Cloudy 
version 08.00 (see the Hazy documentation of Cloudy for 
details) . 

However, since our model can track the abundances of indi- 
vidual elements we know the abundances in each cooling 
halo. In principle, we could recompute a cooling/heating 
rate for each halo using its specific abundances as input 
into Cloudy. This is computationally impractical however. 
Instead, we follow the approach of Martmcz-Sorr ano et al.| 
(2008 1 who perform a principal components analysis (PCA) 
to find the optimal linear combination of abundances which 
minimizes the variance between cooling/heating rates com- 
puted using that linear combination as a parameter and a 
full calculation using all abundances. The best linear combi- 
nation turns out to be a function of temperature. We there- 
fore track this linear combination of abundances at 10 dif- 
ferent temperatures for all of the gas in our models and 
use it instead of metallicity when computing cooling/heating 
rates. 



Compton Cooling: Cole et al. (2000 1 allowed hot halo 



gas to cool via two-body coUisional radiative processes. How- 
ever, as we go to higher redshifts the effect of Compton cool- 
ing must be considered. The Compton cooling timescale is 
given by ( |Peebles 1968) : 



'^Compton 



3mec(l + l/So) 



(31) 



13 As defined by |CoIe et"ar] | |2000| their eqn. Cll) and equal to 
the product of the bulge half-mass radius and the circular velocity 
at that radius. 



where Xc = Wc/nt, is the electron number density, nt 
is the number density of all atoms and ions, Tcmb is the 
cosmic microwave background (CMB) temperature and Tc 
is the electron temperature of the gas. 

The electron fraction, Xc, is determined from photoion- 
ization equilibrium computed using Cloudy (see above). 

Molecular Hydrogen Cooling: The molecular hydrogen 
cooling timescale is found by first estimating the abundance, 
fH2,c, of molecular hydrogen that would be present if there 
is no background of i^2-dissociating radiation from stars. 
For gas with hydrogen number density nn and temperature 



Tv the fraction is (Tegmark et al. 19971 

/h2.c 



3.5 X 10"''r3 '^^[l + (7.4 X 10^(1 + zf-'^^ 



xexp {-3173/(1 + z)}/nHi) 



(32) 



where T3 is the temperature Tv in units of lOOOK and nni 
is the hydrogen density in units of cm"^. Using this initial 
abundance we calculate the final H2 abundance, still in the 
absence of a photodissociating background, as 



/h2 = /Hs.cCXP 



(- 



(33) 



V 51920 A' 

where the exponential cut-off is included to account for col 



lisional dissociation of H2, as in Benson et al. (20061. 

Finally, the cooling time-scale due to molecular hydro- 



gen was computed using ( Galli & Palla 1998 1 

TU, = 6.56419-3'Te/H>Hl'AH2' 

where 

Alte 



Ah2 = 
where 

and 



1 + n" /nn ' 



Ah2(LTE) 
Ah2 [nn 0] ' 



(34) 



(35) 



(36) 



logio Ah2 [nn ^ 0] = -103 + 97.59 ln(T) - 48.05 ln(r)^ 
+10.8 \n{Tf - 0.9032 In(r)'' (37) 

is the cooling function in the low density limit (independent 



of hydrogen density) and we have used the fit given by Galli 
fc Palla (1998}, 



At 



= Ar -I- A„ 



(38) 



is the cooling function in local thermodynamic equilibrium 
and 



Ar = 



A„ 



1 

+3 X lO'^^exp 



■ exp 



0.13 



nui I 



6.7 X 10"^%xp 



( 



3 -1 

ergs cm s , 



5.86 



+1.6 X 10" 



exp 



11.7 



)} 



3 -1 
ergs cm s 



(39) 



(40) 



are the cooling functions for rotational and vibrational tran- 



sitions in H2 ( Hollenbach & McKee 1979 \ 



The model also allows for an estimate of the rate of 
molecular hydrogen formation on dust grains using the ap- 



proach of Cazaux & Spaans (20041. In this case we have to 
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modify equation (13) of Tegmark et al. (19971, which gives 



the rate of change of the H2 fraction, to account for the dust 
grain growth path. The molecular hydrogen fraction growth 
rate becomes: 



/ = fcd/(l -X- 2/) + kmn{l - X - 2/)x, 



(41) 



where / is the fraction of H2 by number, x is the ionization 
fraction of H which has total number density n, 



fed = 3.025 X 10" 



0.01 



lOOK 



(42) 



is the dust formation rate coefficient iCazaux & Spaans 
2004 eqn. 4), and fcm is the effective rate coefficient for H2 



formation (Tegmark et al. 1997 eqn. 14). We adopt the ex- 
pression given by Cazaux & Spaans ( 2004 eqn. 3) for the H 
sticking coefficient, S'h(T') and = 0.53Z for the dust-to- 
gas mass ratio as suggested by Cazaux fc Spaans (2004| ) and 
which results in ^d ~ 0.01 for Solar metallicity. This equa- 
tion must be solved simultaneously with the recombination 
equation governing the ionized fraction x. The solution, as- 
suming x{t) = Xo/{l + xonkit) and 1 — x — 2/ ~ 1 as do 
Tegmark et al. (1997| ), is 



fit) = /o 



exp 



Td 



Ei 



Tr+t 
Td 



)} 



(43) 



where = l/xo/nu/ki, Td = l/nn/fed, ki is the hydrogen 
recombination coefficient and Ei is the exponential integral. 

2.7 Sizes and Adiabatic Contraction 

The angular momentum content of galactic components is 
tracked within our model, allowing us to compute sizes for 
disks and bulges. We follow the same basic methodology as 
Cole et al. (2000 1 — simultaneously solving for the equilib- 



rium radii of disks and bulges under the infiuence of the 
gravity of the dark matter halo and their own self-gravity 
and including the effects of adiabatic contraction — but treat 
adiabatic contraction using updated methods. 

For the bulge component with pseudo-specific angular 
momentum the half-mass radius, rb, must satisfy 



jb = G[Mh(rb) + Md(rb) + Mb(rb)]rb, 



(44) 



where Mh(r), Md(r) and Mb(r) are the masses of dark mat- 
ter, disk and bulge within radius r respectively, and which 
we can write as 



Cb 



[A/h(r-b) + Md(rb) + Mb(rb)]rb, 



(45) 



where Cb = Jb/G. In the original Blumenthal et al. (19861 



treatment of adiabatic contraction the right-hand side of 



eqn. ( 45 1 is an adiabatically conserved quantity allowing us 



to write 



Cb 



M^(rb,o)rb,o, 



(46) 



where is the unperturbed dark matter mass profile and 
rb,o the original radius in that profile. This allows us to 
trivially solve for rb,o and M^{rh,o) and so, assuming no shell 
crossing, Mh(rb) = fh.M^{rh,o), where /h is the fraction of 
mass that remains distributed like the halo. Given a disk 
mass and radius this allows us to solve for rb. 

In the [Gnedin et al. (2004[ ) treatment of adiabatic con- 
traction however, Af(r)r is no longer a conserved quantity. 



Instead, M{{r))r is the conserved quantity where {T)/r-^ 
^ac(r/rh)"'''=. In this case, we write 



= ('■(,) = ^acrh(rb/rh)'' 



Equation ( 45 1 then becomes 



= [A4((r;)) + A/d({r;)) + hH{r'^))]rl 
where 



Cb = 



Cb 



(47) 
(48) 

(49) 



The right-hand side of eqn. ( |48[ ) is now an adiabatically con- 
served quantity and we can write 



c; = A4°(«,o»^b,o. 



(50) 



If we know cj, this expression allows us to solve for r^.o and 
MiAKfi)) which in turns gives Mh(rb) = hM^AKfi))- Of 
course, to find cj, we need to know rb. This equation must 
therefore be solved iteratively. In practice, for a galaxy con- 
taining a disk and bulge, the coupled disk and bulge equa- 
tions must be solved iteratively in any case, so this does not 
significantly increase computational demand. 
The disk is handled similarly. We have 



j1 



= G 



A/h(rd) + yMd-HMb(rd; 



(51) 



where feh gives the contribution to the rotation curve in the 
mid-plane and kd relates the total angular momentum of 
the disk to the specific angular momentum at the half-mass 



radius ( Cole et al. 2000j ). This becomes 
Cd,2 = Afb°((^»r-d,o, 



where 



Cd,2 = 



Cd,2 
^ac 



and 



Cd,2 



Gkl 



(y - ' 1 rdMd. 



(52) 



(53) 



(54) 



This system of equations must be solved simultaneously to 
find the radii of disk and bulge in a given galaxy. Once these 
are determined, the rotation curve and dark matter density 
as a function of radius are trivially found from the known 
baryonic distribution, pre-coUapse dark matter density pro- 
file and the adiabatic invariance of M((r))r. 



2.8 Substructures and Merging 

N-body simulations of dark matter halos have convincingly 
shown that substructure persists within dark matter ha- 



los for cosmological timescales (Moore et al. 19991. More- 



over, recent ultra-high resolution simulations (Kuhlen et al. 
|200£ Bpringcl ct al. 200£ Stadel et al. 2009) > demonstrate that 
multiple levels of substructure (e.g. sub-sub-halos) can ex- 
ist. This "substructure hierarchy" is often neglected in semi- 
analytic models when merging is being considered. For ex- 



ample, Cole et al. (2000 1 and all other semi-analytic models 



14 Andrew J. Benson & Richard Bower 



to dat consider only one level of substructure — a sub- 
structure in a group halo which merges into a cluster imme- 
diately becomes a substructure of the cluster for the pur- 
poses of merging calculations. This is unrealistic and may; 

(i) neglect mergers between galaxies in substructures 
which Angulo et al. (2009) have recently shown to be im- 



portant for lower mass subhalos; 

(ii) bias the estimation of merging timescales for halos 
(and their galaxies). 



Angulo et al. (2009 1 examine rates of subhalo-subhalo merg- 
ers in the Millennium Simulation and find that for subhalos 
with masses below 0.1% the mass of the main halo mergers 
with other subhalos become equally likely as a merger with 
the central galaxy of the halo. They also find that subhalo- 
subhalo mergers tend to occur between subhalos that were 
physically associated before falling into the larger potential. 
This suggests that a treatment of subhalo-subhalo mergers 
must consider the interactions between subhalos and not 
simply consider random encounters as was done, for exam- 



ple, by Somerville & Primack (1999 1 



We therefore implement a method to handle an arbi- 
trarily deep hierarchy of substructure. We refer to isolated 
halos as substructures (i.e. not substructures at all), sub- 
structures of S" halos are called substructures and sub- 
structures of S" halos are 5*"^^ substructures. When a halo 
forms it is an substructure, and when it first becomes a 
satellite it becomes a substructure. 

For 5" substructures with n > 2 we check at the end 
of each timestep whether the substructure has been tidally 
stripped out of its 5*""^ host. If it has, it is promoted to 
being a 5*""^ substructure in the S"~^ substructure which 
hosts its S"~^ host. 



parameters of each subhalo we assume that it can be well 
described as an adiabatic procesf^ As such, the azimuthal 
and radial actions of the orbits: 

I / 



2tv 







and 



Jr = 



rdr, 



(55) 



(56) 



should be conserved (assuming a spherically symmetric po- 
tential). Therefore, at each timestep, we compute Ja and Jr 
for each satellite from the known orbital parameters in the 
current host halo potential. We assume these quantities are 
the same in the new host halo potential and convert them 
back into new orbital parameters rc{E) and e. 

2.8.3 Ttdal Stripping of Dark Matter Substructures 

Given orbital parameters, rc{E) and e we can compute the 
apocentric and pericentric distances of the orbit of each sub- 
halo. At the end of each timestep, for each subhalo we find 
the pericentric distance and compute the tidal field of its 
host halo at that point: 



d 

drh 



GMh(rh) 



(57) 



where uj is the orbital frequency of the subhalo, and find the 
radius, rg, in the subhalo at which this equals 



GM,{r,) 



(58) 



This gives the tidal radius, r^, in the subhalo. 



2.8.1 Orbital Parameters 

When a halo first becomes an subhalo it is assigned 
orbital parameters drawn from the distribution of Benson 



1(200 5) which was measured from N-body simulations. This 
distribution gives the radial and tangential velocity compo- 
nents of the orbit. For later convenience, we compute from 
these velocities the radius of a circular orbit with the same 
energy as the actual orbit, rc(-E), and the circularity (the 
angular momentum of the actual orbit in units of the angu- 
lar momentum of that circular orbit), e. These are computed 
using the gravitational potential of the host halo. 

2.8.2 Adiabatic Evolution of Host Potential 

As a subhalo orbits inside of a host halo the gravitational po- 
tential of that host halo will evolve due to continued cosmo- 
logical infall. To model how this evolution affects the orbital 



^* [Taylor fc Babul (2004[ |, who describe a model of the orbital 
dynamics of subhalos, do account for the orbital grouping of sub- 
halos arriving as part of a pre-existing bound system (i.e. when a 
halo becomes a subhalo its own subhalos are given similar orbits in 
the new host). However, as noted by [Taylor &: Babul (2005l l, they 
do not include the self-gravity of subhalos and so sub-subhalos 
do not remain gravitationally bound to their subhalo. As such, 
sub-subhalos will gradually disperse and cannot merge with each 
other via dynamical friction. 



2.8.4 Promotion through the hierarchy 

After computing tidal radii, for each S-^ subhalo we com- 
pute the apocentric distance of its orbit and ask if this ex- 
ceeds the tidal radius of its host. If it does, the subhalo is 
assumed to be tidally stripped from its host halo and pro- 
moted to an orbit in the host of its host: S" — >■ S"~^ . To 
compute orbital parameters of the satellite in this new halo 
we determine its radius and velocity at the point where it 
crosses the tidal radius of its old host. These are added vec- 
torially (assuming random orientations) to the position and 
velocity of its old host at pericentre in the new host. From 
this new position and velocity values of rc{E) and e are 
computed. 

This approach can handle an arbitrarily deep hierarchy 
of substructure. In practice, the actual depth of the hierar- 
chy will depend on both the mass resolution of the merger 
trees used and the efficiency of tidal forces to promote sub- 
structures through the hierarchy. Given the resolution of the 
trees used in our calculations we find that most substruc- 
tures belonge to the and levels. However, the deepest 
substructure level that we have found at 2; = is S^. 

Halos are expected to grow on the Hubble time, while the 
characteristic orbital time is shorter than this by a factor of \/A 
where A is the overdensity of dark matter halos. This expected 
validity of the adiabatic approximation has been confirmed in 



N-body simulations by Book et al. (2010 1 
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2.8.5 Dynamical Friction 



We adopt the fitting formula found by Jiang et al. (2008 1 



to estimate merging timescales for dark matter substruc- 
tures (and, consequently, the galaxies that they contain). 
The multiple levels of substructure hierarchy in our model 
allow for the possibility of satellite-satellite mergers. We in- 
tend to compare results from our model with N-body mea- 
sures of this process in a future work. 

When a halo first becomes a satellite, we set a di- 
mensionless merger clock, xuf = 0. On each subsequent 
timestep, xuf is incremented by an amount At/rDP where 
tdf is the dynamical friction timescale for the satellite in 
the current host halo according to the expression of [Jiang] 
et al. (2008 1, including the dependence on rc{E). When 



sdf = 1 the satellite is deemed to have merged with the 
central galaxy in the host halo. 

When a satellite is tidally stripped out of its current or- 
bital host and promoted to the host above it in the hierarchy 
the merging clock is reset so that dynamical friction calcula- 
tions start anew in this new orbital host. This is something 
of an approximation since the dynamical friction timescale 
of Jiang et al. (2008] ) is calibrated using satellites that enter 
their halo at the virial radius. As such, it does not explore 
as a sufficiently wide range in rc as is required for our mod- 
els. Furthermore, when promoted to a new orbital host, a 
satellite will have already lost some mass due to tidal effects. 
This is not accounted for when computing a new dynamical 
friction timescale however, and so may cause us to underes- 
timate merging timescales somewhat. 

Dynamical friction also affects the orbital parameters of 
each subhalo. To simplify matters we follow [Lacey fc Cole| 
'{T993I 1 and examme the evolution of these quantities in an 
isothermal dark matter halo. In such a halo, and for a cir- 
cular orbit, rc evolves as 



rc 
rc,o 



= 1 



t 

Tdf 



Therefore, after each timestep we update 
At 



2 
rc 



2 
rc 



2 



(59) 



(60) 



The fractional change in e is assumed to be given by 
{i/e)/{rc/rc) as computed for the current orbit using the 



expressions of Lacey & Cole (19931. This is a function of e 
only and is plotted in Fig. 4 Note that the timescale, tdf, 



used here is that from Jiang et al. (2008| ) and not the one 
from [Lacey fc Cole {1993f. 



2.9 Ram pressure and tidal stripping 



We follow Font et al. (2008 \ and estimate the extent to which 
ram pressure from the hot atmosphere of a halo may strip 
away the hot atmosphere of an orbiting subhalo. In addition, 
we also consider tidal stripping of this hot gas and both ram 
pressure and tidal stripping of material from galaxies. 

Ram pressure and tidal forces are computed at the peri- 
centre of each subhalo's orbit, which we n ow c ompute self- 
consistently with our orbital model (see S2.8l. For an S", 
where i > 1, subhalo we compute the ram pressure force 
from all halos higher in the hierarchy and take the max- 
imum of these to be the ram pressure force actually felt. 




Figure 4. The ratio (i/^)/{rc/rc) for isothermal halos. This 
ratio is used in solving for the evolution of orbital circularity 
and orbital radius under the influence of dynamical friction as 
described in §2.8.5| 



The tidal field (i.e. the gradient in the gravitational force 
across the satellite) includes the centrifugal contribution at 
the orbital pericentre and is given by: 



T 



d GA/(<r) 



AR r2 
The ram pressure is taken to be 

-Pram ~ Phot , host Vorbit 



(61) 



(62) 



where phot, host is the density of hot gas in the host halo at 
the pericentre of the orbit and l^rbit is the orbital velocity 
of the satellite at that position. 



2. 9. 1 Stnpptng of hot halo gas 

We find the ram pressure radius in the hot halo gas by solv- 
ing 

GMeat(rr) 



-Pram — Q^r 



-phot,sat(7^r) 



(63) 



for rr, where Oram is a parameter that we set equal to 2 



as suggested by McCarthy et al. (20081. Similarly, a tidal 



radius is found by solving 



_ 3 G-Msat 
— Q'tidal 5" 



(64) 



for rt , where Qtidai is a parameter that we set equal to unity. 
Once the minimum of the ram pressure and tidal stripping 



radii has been determined we follow [Font et al. (2008"| ) and 
compute the cooling rate of the remaining, unstripped gas 
by cooling only the gas within the stripping radius and as- 
suming that stripping does not alter the mean density of gas 
within this radius. We implement this by giving the satel- 
lite a nominal hot gas mass M^^^ — Mhot + -Mstrip (where 
Mhot is the true hot gas content of the halo) and applying 
the same cooling algorithm as that used for central galaxies 
(except limiting the maximum cooling radius to rstrip rather 
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than -Rv)- This step ensures self-consistency in the treatment 
of the gas cooHng between stripped and unstripped galax- 
ies, and therefore that the colours of satellites are predicted 
correctly. 

The initial stripping of re-heated gas is the same as for 
the hot gas, i.e. the same fraction is transferred from the re- 
heated gas of the satellite to the re-heated gas reservoir of 
the parent halo. We follow Font ct al. (2008 ) in modelling the 
time-dependence of the hot gas mass in the satellite halo and 
refer the reader to that paper for full details. This process 
introduces one free parameter, tstrip which represents the 
time averaged stripping rate after the initial pericentre. We 
treat eatrip as a free parameter which we will adjust to match 
observational constraints. 

The stripping of satellites is also affected by the growth 
of the halo in which the satellite is orbiting. |Font et al.| 
|(2008[ | took this effect into account by assigning each satel- 
lite galaxy new orbital parameters and deriving a new strip- 
ping factor every time the halo doubles in mass compared to 
the initial stripping event. In the present work we directly 
follow the evolution of the pericentric radius and velocity of 
each satellite due to both dynamical friction and host halo 
mass growth. For this reason, we take a different approach 



from Font et al. (20081, computing a new ram pressure ra- 



dius in each timestep instead of only at every mass doubling 
event. 

Any material stripped away from the subhalo is added 
to the halo which provided the greatest ram pressure force. 
For tidal forces, we consider only the contribution from the 
current orbital host as typically if this were exceeded by 
the tidal force from a parent higher up in the hierarchy the 
subhalo would have already been tidally stripped from this 
orbital host and promoted to a higher level in the hierarchy. 



2.9.2 Stripping of galactic gas and stars 

The effective gravitational pressure that resists the ram pres- 



sure force in the disk plane is (for an exponential disk; Abadi 



et al. 19991: 



p — 

erav — 



GMdMg 



'"(I) -'(I) -'■(!) 'Mi) 



,(65) 



where x = r/r^ and Jo, Ii, Kq and Ki are Bessel functions. 
The ram pressure radius is found by solving for the radius 
at which Pgrav = Pthui, where Pram is given by eq. (62 1. 



We assume that any stars in the galaxy which lie beyond 
the computed tidal radius and any gas which lies beyond 
the smaller of the tidal and ram pressure radii are instan- 
taneously removed. Stars become part of the diffuse light 
component of the halo (i.e. that which is known as intra- 



cluster light in clusters of galaxies; see S 4.12.2 1, while gas 
is added to the reheated reservoir of the host halo. The re- 
maining mass of each component (cold gas, disk and bulge 
stars) is computed and the specific angular momentum of 
the remaining material is computed assuming a flat rota- 
tion curve: 



Jdisk 



JdiskO 



E{R)RMR 




E 0.01 r 



40 60 
Step Number 



Figure 5. The remaining mass fraction in an exponential disk 
in a potential giving a flat rotation curve (and ignoring the disk 
self-gravity) subjected to tidal truncation at radius rt/r^ q = 0.1, 
0.3, 1.0, 3.0 and 10.0 (from lower to upper lines) after a given 
number of steps according to our model. The remaining mass 
fraction quickly converges to a near constant value. 
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(66) 



X {Ml - (1 + x,)e-^*] + /g[l - (1 + x^)e-=^-p7) 

for the disk (the last line assuming an exponential disk) 
where R^, = nidai, -Rg = min(rtidai, rram), a;* = R*/Rd, 
xg = i?g/i?d, /. = M./{M, + Ms) and /g = MJ{Nh + M^), 
and 



Jsph JsphO 



/;--^°' p.{R)R^dR/ p.{R)R^dR 

/;"^"' P.{R)RMR/ p.{R)RMR 



(68) 



for the bulge (and which must be evaluated numerically). 
Here, jdisko and jspho are the pre-stripping specific angu- 
lar momenta of disk and spheroid respectively, E*(_R) and 
Egas are the surface density profiles of stars and gas in the 
disk prior to stripping and p*{R) is the stellar density pro- 
file in the spheroid prior to stripping. Since Galform al- 
ways assumes a de Vaucouler's spheroid and an exponen- 
tial disk with stars tracing gas the stripped components will 
readjust to these configurations with their new masses and 
angular momenta. This is, therefore, an approximate treat- 
ment of stripping. In particular, some material will always 
"leak" back out beyond the stripping radius and so is easily 
stripped on the next timestep. Figure [5] demonstrates that 
this is not a severe problem, with the remaining mass frac- 
tion asymptoting to a near constant value after just a few 
steps. 
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2.10 IGM Interaction 



Benson et al. (2002b I introduced metiiods to simultaneously 
compute the evolution of the IGM and the galaxy population 
in a self-consistent manner such that emission from galaxies 
ionized and heated the IGM which in turn lead to suppres- 
sion of future galaocy formation. A major practical limitation 
of [Benson et al.[ s ( |2002b| > method was that it required Gal- 
FORM to be run to generate an emissivity history for the 
Universe which was then fed into a model for the IGM evo- 
lution. The IGM evolution was used to predict the effects 
on galaxy formation and Galform run again. This loop was 
iterated around several times to find a converged solution. 
This problem was inherent in the implementation due to the 
fact that Galform was designed to evolve a single merger 
tree to z = then move onto the next one. 

To circumvent this problem, we have adapted Galform 
to allow for multiple merger trees to be evolved simultane- 
ously: Each tree is evolved for a single timestep after which 
the IGM evolution for that same timestep is computed. This 
allows simultaneous, self-consistent evolution of the IGM 
and galaxies without the need for iteration. 

The model we adopt for the IGM evolution is essen- 
tially identical to that of [Benson et al. (2002b[ ), and consists 
of a uniform IGM (with a clumping factor to account for en- 
hanced recombination and cooling due to inhomogeneities) 
composed of hydrogen and helium and a photon background 
supplied by galaxies and AGN. The reader is therefore re- 



ferred to Benson et al. (2002b I for a full discussion. Here we 
will discuss only those aspects that are new or updated. 



2.10.1 Emissivity 

The two sources of photons in our model are quasars and 
galaxies. For AGN we assume that the spectral energy dis- 



tribution (SED) has the following shape (Haardt & Madau 



19961 



/.(A) « 




if A < 1216A; 

if 1216A < A < 2500A; 

if A > 2500A, 



(69) 



where the normalization of each segment is chosen to give 
a continuous function and unit energy when integrated over 
all wavelengths. The emissivity per unit volume from AGN 
is then 



EAGN = /csc,AGNe.p.C^/,^(A), 



(70) 



where e, =0.1 is an assumed radiative efhciency for accre- 
tion onto black holes, p. is the rate of black hole mass growth 
per unit volume computed by Galform and /osc.agn is an 
assumed escape fraction for AGNphotons which we fix at 
10~^ to produce a reasonable epoch of Hell reionization. 

The emissivity from galaxies was calculated directly by 
integrating the star formation rate per unit volume predicted 
by Galform over time and metallicity to give 



tgal 



/cBc.gal(t')Af*(t', ^)i.(tnow ^ t' , Z[t'])At' , (71) 



where M*(f, Z) is the rate of star formation at metallicity Z, 
Lv{t, Z) is the integrated luminosity per unit frequency and 
per Solar mass of stars formed of a single stellar population 



of age t and metallicity Z and /csc.gai is the escape fraction 
of ionizing photons from the galaxy. 

The fraction of ionizing photons able to escape from 
the disk of each galaxy is computed using the expressions de- 



rived by Benson et al. (2002a I (their eqn. A4) which is a gen- 
eralization of the model of [Dove & ShuU (1994 1 in which OB 
associations with a distribution of luminosities ionize holes 
through the neutral hydrogen distribution through which 
their photons can escape. 

The sum of eagn and egai gives the number of photons 
emitted from the galaxies and quasars in the model. 



2.10.2 IGM Ionization State 

The ionization state of the IGM is computed just as in 



Benson et al. (2002b I except that we use effective photo- 



ionization cross-sections that account for the effects of sec- 
ondary ionizations and a re given by [ShuU fc van Steenberg| 
( 1985 as re-expressed by Venkatesan et al. 2001 1: 



^'n{E) = 1 



^c(£) 



E-E^ 



E- E 



+ ( 1 + <^Hcl 

E 



Eh 

E — Euc 



1 



Euc 

Ehc 



"l9.95eVy 

crHo(-B) 



au{E) 



+ 



-Ehc 
E — Eh 
' 24.6 



o-Hc(£) 



aniE) 



(72) 



(73) 



where (j{E) is the actual cross section (Verner & Yakovlev 
p95| ) and 

().4()92n1.7592 



</>Hi = 0.3908(1 
fnc = 0.0246(1 - X. 
feci = 0.0554(1 - X. 



0.4049n1.6594 



0.4614n1.6660 



(74) 
(75) 
(76) 



2.10.3 IGM Thermal State 



Heating of the IGM is treated as in Benson et al. (2002b \ 
with the exception that we account for heating by secondary 
electrons. Photoionization heats the IGM at a rate of 



■Jphoto 



{E - Ei)ca' (E)nin-,{E)£dE 



(77) 



where Ei is the energy of the sampled photons which is asso- 
ciated with atom/ion number density n^, c is speed of light, 
a' is the effective partial photo-ionization cross section (ac- 
counting for secondary ionizations) for the ionization stages 
of H and He, nji^E) is the number density of photons of 
energy E, Ei is the ionization potential of i and index i rep- 
resent the different atoms and ions, H, H"*", He, He"^ and 
He^"*". In the above, £ accounts for heating by secondary 



electrons and is given by ( ShuU & van Steenberg 1985 1 
£ = 0.9971[1 - (1 - a:" '''^^)'-^'''^]. 



(78) 



2.10.4 Suppression of Baryonic Infall into Halos 

According to [Okamoto et al. (2008[ ), the mass of baryons 
which accrete from the IGM into a halo after reionization is 
given by 



Mi + Ma, 



(79) 
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where 



= 2^ exp 



St 



(80) 



and where the sum is taken over the progenitor halos of the 
current halo, St is the time since the previous timestep and 
tevp is the timescale for gas to evaporate from the progenitor 
halo and is given by 



to 



Rh/c^A 

cvp J 

OO 



if 7v *C Tcvp ; 

if Tv ^ Tcvp ■ 



(81) 



Here, Tevp is the temperature below which gas will be heated 
and evaporated from the halo. We follow [Okamoto et al.| 
(2008 1 and compute Tevp by finding the equilibrium temper- 
ature of gas at an overdensity of Aevp = 10'^. The accreted 
mass Afacc is given by 



no 




if Ti^ir ^ Ta.c,c 
if Tvir Tacc 



(82) 



where Tacc is the larger of the temperature of IGM gas adi- 
abatically compressed to the density of accreting gas and 
the equilibrium temperature, Tcq, at which radiative cool- 
ing balances photoheating for gas at the density expected at 
the virial radius. This ensures that a sensible temperature is 
used even when the photoionizing background is essentially 
zero. 

The value of Tacc is computed at each timestep by 



searching for where the cooling function (see S 2.6.7 1 crosses 
zero for the density of gas just accreting at the virial radius 
(for which we use one third of the halo overdensity; |Okamoto| 



et al. 20081. 



2.11 Recycling and Chemical Evolution 



In Cole et al. (2000] ) the instantaneous recycling approxi- 
mation for chemical enrichment was used. While this is a 
reasonable approximation for z = it fails for high redshifts 
(where the main sequence lifetimes of the stars which do 
the majority of the enrichment become comparable to the 
age of the Universe). It also prevents predictions for abun- 
dance ratios (e.g. [a/Fe]) from being made and ignores any 
metallicity dependence in the yield. 



Nagashima et al.l (|2005a see also Nagashima et al 



2005b Arrigoni et al. 2009 1 previously implemented a non- 



instantaneous recycling calculation in Galform. We imple- 
ment a similar model here, following their general approach, 
but with some specific differences. 

The fraction of material returned to the ISM by a stellar 
population as a function of time is given by 



R{t) 



[M - Mr(M;Z)](t>(M) 



M{t;Z) 



AM 



(83) 



where 4>{M) is the initial mass function (IMF) normalized 
to unit stellar mass, Mr{M) is the remnant mass of a star of 
initial mass M. Here, M{t) is the mass of a star with lifetime 
t. Similarly, the yield of element i is given by 



M,{Mo;Z)4,{Mo 



M{t;Z) 



dMo 



(84) 



where Mi{Mo; Z) is the mass of metals produced by stars of 
initial mass AIq. For a specified IMF we compute R{t; Z) and 



yi{t; Z) for all times and elements of interest. This means 
that, unlike most previous implementations of Galform, 
the recycled fraction and yield are not free parameters of 
the model, but are fixed once an IMF is chosen. However, 
it should be noted that significant uncertainties remain in 
calculations of stellar yields, which may therefore influence 



our calculations. Note that, unlike Nagashima et al. (2005a I, 



we include the full metallicity dependence in these functions. 



Stellar data are taken from Portinari et al. (1998) for low and 
intermediate mass stars and Marigo (2001 1 for high mass 
stars. 

In Galform the evolution of gas and stellar masses in 
a galaxy are controlled by the following equations'^ 



Mr 



Me, 



= -(! + /?') 



+ Mr + Mn 



where 



Tdisk 



/"dyn^bulg 



(200km S 



for disks 
for bursts. 



(85) 
(86) 

(87) 



is the star formation timescale, Tdisk is the dynamical time 
at the disk half-mass radius, rbuigc is the dynamical time at 
the bulge half-mass radius, /dyn — 2 and f5' quantifies the 
strength of supernova feedback (see ^2.12[ ). In [Cole et al.| 
(20001, the instantaneous recycling approximation implies 
that Mr oc Mgas/i"*, and the cosmological infall term Mintaii 
is approximated as being constant over each short timestep. 
This permits a simple solution to these equations. In our 
case, we retain the assumption of constant Minfaii and fur- 
ther assume that the mass recycling rate, Mr, can be ap- 
proximated as being constant throughout the timeste]J^ 
We therefore write 



Mr = 



A/R,past + Mr,, 

At 



where At is the timestep, 

pto+At 

MR,past = / dt" 
J to 



dt'M,{t')R{t" -t') 



(89) 



is the mass of gas returned to the ISM from populations of 
stars formed in previous timesteps (and is trivially computed 
from the known star formation rate of the galaxy on past 
timesteps) and 



Mr, now — 



dt'' 



At'M,{t')R{t" -t'), 



(90) 



is the mass returned to the ISM by star formation during 
the current timestep. With these approximations, the gas 
equations always have the solution 



Mgas(t) = MgasO exp ) + A/input rcft■ 
-^cfi■ 



1 — exp ( — 



Tcff 



^® These are identical to those given in |Cole et al.| |2000[ their 
equations 4.6 and 4.8) except for the explicit inclusion of the 
recycling terms — [Cole et al. (20do| | included these using the in- 
stantaneous recycling approximation. 

This will be approximately true if the timestep is sufficiently 
short that RM < R. 
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where Mgaso is the mass of gas at time t 
the start of the timestep and 



(measured from 



Minfall 



Ml 



R ,past 



Tcff 



At 



/Rl(At, Tcff) 



Mr, past 
At 



-/Ro(At)| 

X {(1 + /?) + [/ri (At, Toff) - lKoiAt)]/At}(±) 



where 



Ino{t) = 



/Ri(t,r) = 



R{t - t')dt' 



exp{-t'/T)R{t-t')dt'. 



(93) 



(94) 



In the above, the effective e-folding timescale for star for- 
mation (accounting for supernovae (SNe) driven outflows), 
Teff, is given by 



Toff = 



1 + 13' 



(95) 



where /3' measures the strength of SNe feedback and is de- 



fined below in eqn. ( 103 1 



The evolution of the metal mass is treated in a similar 
way, assuming a constant rate of input of metals from infall, 
star formation from previous timesteps and star formation 
from the current timestep. Metals in the cold gas reservoir 
of a galaxy are assumed to be uniformly mixed into the gas, 
such that the reservoir has a uniform metallicty. Metals then 
flow from the cold gas reservoir into the stellar phase and out 
into the reheated reservoir at a rate proportional to the star 
formation rate and mass outflow rate respectively, with the 
constant of proportionality being the cold gas metallicity. 
Material recycled from stars to the cold phase carries with 
it metals corresponding to the original metallicity of those 
stars, augmented by the appropriate metal yield. Finally, gas 
infalling from the surrounding halo may have been enriched 
in metals by previous galaxy formation and so deposits met- 
als into the cold phase gas at a rate proportional to the mass 
infall rate, with proportionality equal to the (assumed uni- 
form) metallicity of the notional proflle gas. Apart from the 
fact that metals from stellar recycling and yields are not 
added instantaneously to the cold reservoir this treatment 
of metals remains identical to that of Cole et al. (2000| . The 
net rate of metal mass input to the cold phase (from both 
cosmological infall and returned from stars) is 



Mz, in 



Mz, infall 



+ - 



^]/Rl(At,reff) 



-hoiAt) 



At[{l + p) + {lKi{At, Teff) - Im){At))At] 



2.11.1 Star Bursts 

In previous implementations of Galform star bursts were 
assumed to have an exponentially declining star formation 
rate. Such a rate results from assuming an instantaneous 
star formation rate of 

Mcold 



M* = 



(99) 



where Ti, is a star formation timescale (fixed throughout the 
duration of the burst), an outflow rate proportional to the 
star formation rate and a rate of recycling given by RM^. 
The resulting differential equations have a solution with an 
exponentially declining star formation rate. 

When the instantaneous recycling approximation is 
dropped the rate of recycling is no longer proportional to the 
star formation rate and the differential equations no longer 
have an exponential solution. We choose to retain the orig- 
inal star formation law (eqn. 99 1 and solve the differential 
equations to determine the star formation rate, outflow rate 
etc. as a function of time in the burst. The resulting set 
of equations have solutions identical to those in §2.11| but 
with zero cosmological infall terms. Recycled material and 



the effects of feedback (see S 2.12 1 are applied to the gas in 
the burst during the lifetime of the burst. Any recycling and 
feedback occurring after the burst is flnished are applied to 
the disk. 



In Cole et al. (2000 1 while bursts were treated as having 



finite duration for the purposes of computing the luminosity 
of their stellar populations at some later time, the change 
in the mass of the galaxy due to the burst occurred instan- 
taneously. We drop this approximation and correctly follow 
the change in mass of each component (gas, stars, outflow) 
during each timestep. 



2.12 Feedback 

Feedback from supernovae is also modifled to account for the 
delay between star formation and supernova. In |Cole et al.| 
(2000 1 the outflow rate due to supernovae feedback was 



where 



galaxy 



(100) 



(101) 



Vhot and ahot are parameters of the model (we allow for two 
different values of Vhot, one for quiescent star formation in 
disks and one for bursts of star formation) and Vgaiaxy is 
the circular velocity at the half-mass radius of the galaxy, 
determines the strength of feedback and is a function of the 
depth of the galaxy's potential well. We modify this to 



-f^]7pi(At, Teff) + ^/po(At) M,^^ = p'M, 



At[{l + + (/pi(At,reff) - Jpo(Ai))At] 

where Mz^R.past is the mass of metal i recycled from star 
formation in previous timesteps and 

ft 



where 



13' ^13 



J^Mt')NsNc{t-t')t' 



(102) 



(103) 



/po(t) = / p{t-t')dt', 
Jo 

/pi(t,r) = / exp{-t'/T)p{t-t')dt' 
Jo 



(97) 
(98) 



where NsNc{t) is the total number of SNe (of all types) aris- 
ing from a single population of stars after time t, such that 
the outflow rate scales in proportion to the current rate of 
supernovae but produces the same net mass ejection after 
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infinite time (for constant In fact, we compute /? using 
the present properties of the galaxy at each timestep. The 
qualifier "(II)" appearing in the quantity Ai'g^^(cxD) in the 



galaxy has come to the forefront (Croton et al. 200( power 



et al. 200qpomerville et al. 2008b[|. We adopt the black hole 
growth model of|Malbon et al. (2007*) and the AGN feed- 



denominator of eqn. ( 103 1 indicates that we normalize the back model of Bower et al. (2006 j as modified by Bower 



outflow rate by reference to the number of supernovae from et al. (2008 1. The reader is referred to those papers for a full 



our adopted Population II IMF (see S 2.14 1. This results in 
the outflow correctly encapsulating any differences in the ef- 
fective number of supernovae between Population II and III 
stars. For supernovae rates, we assume that all stars with 
initial masses greater than 8Mq will result in a Type II su- 
pernova allowing the rate to be found from the lifetimes of 
these stars and the adopted IMF. We adopt the calculations 



of [Nagashima et al. (2005a[ ) to compute the Type la SNe spectral synthesis library 



description of our implementation of AGN feedback. 



2.14 Stellar Populations 

We consider both Pop II and Pop III stars. To compute 
luminosities of Population II stellar populations we employ 
the most recent versioi p of the Conroy, Gu nn & White 



rate. 



Since /3' appears in the gas equations of j 2.11 but 
also depends on the star formation rate during the current 
timestep we must iteratively seek a solution for /3' which is 
self-consistent with the star formation rate. We find that a 
simple iterative procedure, with an initial guess of /?' = /3 
quickly converges. 

When gas is driven out of a galaxy in this way it can 
be either reincorporated into the Mreheatcd reservoir in the 
notional hot gas profile of the current halo, or it can be 
expelled from the halo altogether and allowed to reaccrete 
only further up the hierarchy once the potential well has 
become deeper. 

We assume that the expelled fraction is given by 



/oxp = exp - 



X 4>V 



(104) 



such that the rate of mass input to the reheated reservoir is 

= (1 - /cxp)/3'M.. (105) 

Here, is a dimensionless parameter relating the depth of 
the potential well to (we set = 1 always), V is the 
circular velocity of the galaxy disk or bulge (for quiescent 
or bursting star formation respectively) and (e) is the mean 
energy per unit mass of the outflowing material. We further 
assume 

(e) = ^AexpolV^ (106) 

where Aexpei is a parameter of order unity relating the energy 
of the outflowing gas to the potential of the host galaxy, and 
will be treated as a free parameter to be constrained from 
observations (we actually allow for Acxpoi to have different 
values for quiescent and bursting star formation; see ij3|. We 
then proceed to the parent halo and allow a fraction 



/acc = exp 



* max 



exp 



vl 



(107) 



to be reaccreted into the hot gas reservoir of the notional 
proflle, where Knax is the maximum of y/ AcxpciV^ and any 
parent halo Vv yet found. We then proceed to the parent's 
parent and repeat the accretion procedure, continuing until 
the base of the tree is reached. In this way, all of the gas will 
be reaccreted if the potential well becomes sufficiently deep. 



2.13 AGN feedback 

In recent years, the possibility that feedback from AGN- 
plays a significant role in shaping the properties of a forming 



a Chabrier IMF ( Chabrier 2003 1 



Conroy et al. 2009| p^ We adopt 



(t)(M) oc 



exp 

M- 



( 1 [logioA//A/cl^ A 

y 2 ) 



forM < IMe (^Qg^ 
for A/ >IMq, 



where Mc = O.O8M0 and a — 0.69 and the two expressions 
are forced to coincide at IMq. Recycled mass fractions, yield 
and supernovae rates are computed self-consistently from 
this IMF as described in pTl] and pT2l and are shown in 
Fig-i 

For Population III stars (which we assume form be- 
low a critical metallicity of Zait ~ lO"*.^©) we adopt IMF 
"A" from [Tumlinson (2006[ ). Spectral energy distributions 
for this IMF as a function of population age were kindly pro- 
vided by J. Tumlinson. Lifetimes for these stars are taken 
from the tabulation given by |TumUnson et al. (2003] ). Recy- 
cled fractions and yields and energies from pair instability 
supernovae are computed using the data given by |Heger fc| 
[Woosley (2002| ). Recycled mass fractions, yield and super- 
novae rates are computed self-consistently from these Pop- 
ulation III stars as shown in Fig. [6] by green lines. 



2.14-1 Extinction by Dust 



Cole et al. (2000 1 introduced a model for dust extinction 
in galaxies which significantly improved upon earlier "slab" 



models. In Cole et al. (20001 the mass of dust is assumed to 
be proportional to the mass and metallicity of the ISM and 
to be mixed homogeneously with the ISM (possibly with a 
different scale height from the stars) and to have properties 
consistent with the extinction law observed in the Milky 
Way. To compute the extinction of any galaxy, a random 
inclination angle is selected and the extinction computed 
using the results of radiative transfer calculations carried 



out by Ferrara et al. (19991 



Following Gonzalez-Perez et al. (20091, we extend this 
modeQby assuming that some fraction, /cioud, of the dust 
is in the form of dense molecular clouds where the stars form 



Specifically, v2.0 downloaded from 
^http: //www, astro .prlncet on. edu/^cconroy/SPS/^ with bug 
fixes up to January 7, 2010. 

For calculations of IGM evolution we do not use the I Con- 1 
|roy et al. (2009| spectra because they assign stars hotter than 
5 X 10«K pure blackbody spectra. This leads to an unrealistically 
large ionizing flux for young, metal rich populations. We there- 
fore instead use the |Bruzual fc Chariot (2003| l spectral synthesis 
library for IGM evolution calculations. 

'^^ An alternative method for rapidly computing dust extinction 
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Figure 6. Upper, middle and lower panels show the recycled frac- 
tion, yield and effective number of supernovae respectively for a 
Chabrier IMF (two metallicities, defined as the mass fraction of 
heavy elements, are shown: 0.0001 as red lines and 0.0501 as blue 
lines) and for metal free Population III stars with type "A" IMF 
from [Tumlinson (2006[ l (green lines). Top-panel: The fraction of 
mass from a single stellar population, born at time t = 0, recy- 
cled to the interstellar medium after time t. Middle-panel: The 
total metal yield from a single stellar population born at time, 
t = 0, after time t is shown by the solid lines. Dotted and dashed 
lines show the yield of oxygen and iron respectively. Lower-panel: 
Cumulative energy input into the interstellar medium, expressed 
as the number of equivalent supernovae, per unit mass of stars 
formed as a function of time. The dotted line indicates the con- 
tribution from stellar winds, the solid line the contribution from 
Type II supernovae and the dashed line the contribution from 
Type la supernovae. 



(see |Baugh et al. 200^pacey et al. 2010[ ). Stars are assumed 
to form in these clouds and to escape on a timescale of Tquics 
(for quiescent star formation in disks) or Tburst (for star for- 
mation in bursts), which is a parameter of the dust model 



(Granato et al. 20001, so these stars spend a significant 
fraction of their lifetime inside the clouds. Since massive, 
short-lived stars dominate the ultraviolet (UV) emission of 
a galaxy this enhances the extinction at short wavelengths. 

To compute emission from dust we assume a far infrared 
opacity of 



^iCA/Ai)-"! 

'^1 (Abroak/Al^ 



-0 



^ (A/Abrcak) 



for A < Abrcak 
for A > A break, 



(109) 



and re-emission within the Galform-|-Grasil frameworks based 



Table 1. Parameters of the dust model used throughout this 
work. The parameters are defined in j ]2.14.1| 



Parameter 


Value 


/"cloud 


0.25 


^ burst 


1.0 


'tquics 


1 Myr 


^burst 


1 Myr 


-^l.disk 


SOfim 


^breakjdisk 


lOOOO^m 


/3l,disk 


2.0 


/32,disk 


2.0 


-^1, burst 


30/xm 


-^brcak, burst 


lOO^tm 


burst 


1.6 


/32, burst 


1.6 



where the opacity normalization at Ai = 30^m is chosen to 
be Ki — 140cm^/g to reproduce the dust opacity model used 
in Grasil, as described in |Silva et al. (1998| ). The dust grain 
model in Grasil is a slightly modified version of that pro- 



posed by praine fcLee (1984) . Both the |Draine fc Lee (1984 1 
and Grasil dust models have been adjusted to fit data on 
dust extinction and emission in the local ISM (with much 
more extensive ISM dust emission data being used by |Silva| 
|et al. 1998[ ). The normalization is set at 30/im because the 
dust opacity in the Draine & Lee (19841 and Grasil models 
is well fit by a power-law longwards of that wavelength, but 
not shortwards. The dust luminosity is then assumed to be 



gas 1 



(110) 



where B^{T) = [2h!/^/c^]/[exp(hi//kr) - 1] is the Planck 
blackbody spectrum and Mz,gas is the mass of metals in gas. 
The dust temperature, T, is chosen such that the bolometric 
dust luminosity equals the luminosity absorbed by dust. 

Values of the parameters used in dust model are given 
m Table [l] and were found by |Gonzalez-Perez et al. (2009[ ) to 
give the best match to the results of the full Grasil model. 

This extended dust model, including diffuse and molec- 
ular cloud dust components, provides a better match to the 
detailed radiative transfer calculation of dust extinction car- 



ried out by the spectrophotometric code Grasil ( Silva et al. 
199^paugh et al. 200^^augh et al. 200^|:.acey et al. 2008 \ 
while being orders of magnitude faster, although it does 
not capture details such as polycyclic aromatic hydrocar- 
bon (PAH) features. 



Fontanot et al. (2009b I have explored similar models 
which aim to reproduce the results of Grasil using simple, 
analytic prescriptions. They found that by fitting the results 
from Grasil they were able to obtain a better match to 
the extinction in galaxies than previous, simplistic models 
of dust extinction had been able to attain. In this respect, 
our conclusions are in agreement with theirs — the model we 
describe here provides a significantly better match to the 
results of the full Grasil model than, for example, the dust 



extinction model described by Cole et al. (2000 1 



on artificial neural networks is described by Almeida et al. (2009 1 



At high redshifts model galaxies often undergo peri- 
ods of near continuous bursting as a result of experiencing 
disk instabilities on each subsequent timestep. This rather 
chaotic period of evolution is not well modelled presently — 
it is treated as a sequence of quiescent gas accretion periods 
punctuated by instability-triggered bursts while in reality we 
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expect it to correspond more closely to a near continuous, 
high star formation rate mode somewhere in between the 
quiescent and bursting behaviour. While our model proba- 
bly estimates the total amount of star formation during this 
period reasonably well (as it is controlled primarily by the 
cosmological infall rate and degree of outflow due to super- 
novae) we suspect that it does a rather poor job of account- 
ing for dust extinction. After each burst the gas (and hence 
dust) content of each galaxy is reduced to zero, resulting 
in no extinction. Our model therefore tends to contain too 
many dust-free galaxies at high redshifts. To counteract this 
effect we force galaxies in this regime to be observed during 
a bursting phase, so that they always experience some dust 
extinction. 

Dust remains one of the most challenging aspects of 
galaxies to model. We will return to aspects of our model 
related to dust (utilizing the more detailed Grasil model) 
in a future work, but note that even this is unlikely to be 
sufficient — what is needed is a better understanding of the 
complicated distribution of dust within galaxies, particularly 
during these early, chaotic phases. 

Indeed the distribution of star formation within galaxies 
at z = 3 to 5 has recently become within reach of observa- 
tional studies ( Stark et al. 20Qi ^Imegreen et al. 200£ |l/ehnert 



et al. 200!: 3winbank et al. 20091. It seems that this aspect 



of the model is indeed supported by observational data. A 
future project will be to compare the internal properties of 
observed galaxies at these redshifts with those predicted by 
the model. 



2.15 Absorption by the IGM 

Where necessary, we model the attenuation of galaxy SEDs 
by neutral hydrogen in the intervening IGM using the model 



of |Meiksin (2006| ). 



3 MODEL SELECTION 

The model described above has numerous free parameters 
which reflect our ignorance of the details of certain physical 
processes or order unity uncertainties in (e.g. geometrical) 
coefficients. To determine suitable values for these parame- 
ters we appeal to a broad range of observational data and 
search the model parameter space to find the best fit model. 

The problem of how to implement the computationally 
challenging problem of fitting a complicated semi-analytic 
model with numerous free parameters to observational data 
has been considered before by [Henriques et al. (2009[ ) and 
[Bower et al. (20"lO L To constrain model parameters in this 
work we use the "Projection Pursuit" method of |Bower et al.| 
|(2010[ ). We give a brief description of that method here and 
refer the reader to Bower et al. (2010| for complete details. 

Running a single set of model parameters, including all 
of the redshifts and wavelengths required for our analysis, is 
a relatively slow process. In particular, running a model with 
self-consistently computed IGM evolution is entirely imprac- 
tical for a parameter space search. We therefore chose to run 
models without a self-consistently computed IGM or pho- 
toionizing background. Even then, each model takes around 
2 hours to run on a fast computer. To mimic the effects of 
a photoionizing background we adopt the "Vcut-2cut" model 



described by Font et al. (2009| and which they show to re- 
produce quite well the results of the self-consistent calcula- 
tion. Briefiy, this model inhibits cooling of gas in halos with 
virial velocities below Vcut at redshifts below Zcut. We then 
include Vcut and Zcut as parameters in our fitting process. 

This approach is not ideal, but is required due to com- 
putational limitations. [Bower et al. (2010| ) show that local 
(i.e. low redshift) properties of the model are not signifi- 
cantly affected by the inclusion of self-consistent reioniza- 
tion (i.e. those data do not constrain Kut or Zcut), and, 
where they are, the "Vcut-Zcut" model provides a reason- 
able approximation Font et al. (2009). In any case, as we 



will discuss below, some manual tuning of parameters is still 
required after the automated search of parameter space is 
completed. This manual search is then conducted using the 
fully self-consistent IGM calculation. 

We envision the problem in terms of a multi- 
dimensional parameter space into which constraints from 
observational data are mapped. Given the large number of 
model parameters and the fact that running a single realiza- 
tion of the model requires a significant amount of computer 
time, we can not perform a simple grid-search of the pa- 
rameter space on a sufficiently fine grid. Instead, we begin 
by specifying plausible ranges for model parameters. The 
ranges considered for each parameter are listed in Table [2]— 
for some parameters we choose to consider the logarithm 
of the parameter as the variable in our parameter space, 
to allow for efficient exploration of several decades of pa- 
rameter value. We scale each model parameter such that it 
varies between and 1 across this allowed range. We then 
generate a set of points in this limited and scaled model 



parameter space using Latin hypercube sampling (McKay 



et al. 19791, thereby ensuring an efficient coverage of the 



parameter space. A model is run for each set of parameters 
and a goodness of fit measure computed. 

The choice of a goodness of fit measure is important and 
non-trivial (see [Bower et al. 2010] ) . We do not expect our 
model to fit all of the constraints in a statistically rigorous 
manner, as the model is clearly approximate. The Bayesian 
approach to this is issue is to assign a prior assessment of 
the reliability of the model to each of the data set compar- 
isons and to define a correlation matrix reflecting the a priori 
connections between datasets. This concept (referred to as 
"model discrepancy" in the statistical literature) is discussed 
in detail for z = luminosity function constraints in ,Bower| 
[et al. (2010| ). However, in the present paper, we needed a 
simpler approach to the problem. We therefore adopted a 
non-Bayesian methodology of simply summing for each 
dataset that we used. This has the advantage of simplicity, 
but clearly there may be more appropriate choices for the 
relative weighting of different data sets: we will explore this 
issue in a future paper. There is little doubt that a better 
measure of goodness of fit could be found. In particular, 
the relative weightings given to each dataset should really 
reflect how well we think the model performs in that par- 
ticular quantity, how accurately we think that we have been 
able to match any observational selection and, inevitably, 
how much we believe the data itself. These are extremely 
thorny issues to which, at present, we do not have a good 
answer. 

Speciflcally, in this work, the goodness of fit measure is 
taken to be 
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Table 2. The allowed ranges for each parameter in our fitting 
parameter space. For some parameters, we choose to use the log- 
arithm of the parameter to allow efficient exploration of several 
decades of parameter value. 



Parameter 


Minimum 


Maximum 


ho 


0.6750 


0.7270 




0.04320 


0.04920 


Ao 


0.7142 


0.7278 


CS 


0.7650 


0.8690 


ris 


0.9320 


0.9880 


Vcut/km s~l 


10.00 


50.00 


Z^cut 


5.000 


13.00 


loglo("cool) 


-1.523 


0.4771 


log]^Q(Q:romovo) 


-1.523 


0.0000 


loglo{acorc) 


-2.000 


-0.5229 


loglo{«*) 


-3.523 


-1.301 


a* 


-4.000 


1.000 


Vhot,disk/km s-i 


100.0 


550.0 


1 /" /I — 1 

Vhot,burst/km s 


100.0 


550.0 


Q^hot 


1.000 


3.700 


l°glo(-''oxpcl,disk) 


-1.523 


1.000 


l°glO (-^oxpcl, burst ) 


-1.523 


1.000 


loglo(««) 


-2.398 


-1.000 


loglo{^») 


-3.000 


-1.000 




-3.000 


-1.523 


loglo{"rohcat) 


-1.523 


0.4771 


logloC/cllip) 


-2.000 


-0.3010 


loglO (/burst) 


-2.000 


-0.3010 


loglO {/gas, burst) 


-1.523 


-0.3010 


-B/Tburst 


0.0000 


1.000 


^ac 


0.7000 


1.000 


tfac 


0.7000 


1.000 


^d,gas 


0.7000 


1.150 


logloC'^strip) 


-2.000 


0.0000 



~2 Xi 



(111) 



where Xi is the usual goodness of fit measure for dataset i, 
Ni is the number of degrees of freedom in that dataset and 
Wi is a weight assigned to each dataset. The sum is taken 
over all datasets shown in !|4]and, additionally, cosmological 
parameters were allowed to vary within the 2a intervals per- 
mitted by the [Dunkley et al. (2009[ ) constraints, and were 
included in the goodness of fit measure using a Gaussian 
prior. When computing for each dataset we estimate the 
error in each datum to be the sum in quadrature of the 
experimental error and any statistical error present in the 
model due to the finite number of Monte Carlo merger tree 
realizations that we are able to carry out. This ensures that 
two models which differ by an amount comparable to the 
random noise in the models have similar values of x^- The 
specific datasets used, along with the weights assigned to 
them (estimated using our best judgement of the reliability 
of each dataset and the Galform's ability to model it) are 
listed in Tabled 

Once a set of models have been run, a principal com- 
ponents analysis is performed on the goodness of fit values 
of those models with x^ values in the lower 10**^ percentile 
of all models to find which linear combinations of parame- 
ters provide the minimum variance in goodness of fit. These 
are the parameter combinations that are most tightly con- 
strained by the observational data. A principal component 



with low variance implies that this particular combination of 
the parameters is tightly constrained if the model is likely 
to produce an acceptable fit. Of course, even if this con- 
straint is satisfied, a good model is not guaranteed; rather 
we can be confident that if it is not satisfied the fit will 
not be goocp] When analysing the acceptable region in this 
way, we also need to bear in mind that the PCA assumes 



that the relationships are linear, whereas Bower et al. (2010 1 
show that the actual acceptable space is curved. This will 
prevent any of the suggested projections being arbitrarily 
thin and limit the accuracy of constraints. Nevertheless, the 
procedure substantially cuts down the volume of parameter 
space where model evaluations need to be run. These linear 
combinations are used to define rotated axes in the param- 
eter space within which we select a new set of points again 
using Latin hypercube sampling. The process is repeated 
until a suitably converged model is founcj^ This process is 
not fast, requiring around 150,000 CPU hourj^ but does 
produce a model which is a good match to the input data. 

Figure [7| demonstrates the efficacy of our method us- 
ing four 2D slices through the multi-dimensional parameter 
space. The colour scale in each panel shows constraints on 
two of the model parameters, while the projections below 
and to the left of the panel indicate the constraints on the 
indicated single parameter. Contours illustrate the relative 
number of model evaluations which were performed at each 
point in the plane. It can be clearly seen that our "Projec- 
tion Pursuit" methodology concentrates model evaluations 
in those regions which are most likely to provide a good 
fit. The nominal best-fit model is indicated by a yellow star 
in each panel. Despite the large number of models run we 
do not believe that this precise point should be considered 
as the "best" model — the dimensionality of the parameter 
space is so large that we do not believe that it has been suf- 
ficiently well mapped to draw this conclusion. Additionally, 
we also need a model discrepancy matrix — without this, we 
can not say whether a model is acceptable (in the sense that 
it should only agree with the data as well as we expect given 
the level of approximation in the model). Without the dis- 
crepancy term, we will tend to overfit the model. Instead, 
we utilize these results to suggest the region of parameter 



This is only strictly true if the relationships between and 
the parameters are approximately linear and unimodal. If there 
exists a separate small island of good values somewhere, our 
PCA-|-Latin Hypercube method might happen to miss the re- 
gion, or it might not exert sufficient pull on the PCA compared 
to the large region and might be subsequently ignored. The ad- 
vantage of the emulator approach used by [Bower et al. (2010[ l is 
that it gives an estimate of the error made by excluding regions 
from further evaluations. 

In practice these calculations were run on distributed comput- 
ing resources (including machines at the ICC in Durham, Tera- 
Grid and Amazon EC2. Each machine was given an initial small 
set of models to run. After running each model, the results were 
transferred back to a central server. Periodically, the server would 
collate all available results, perform the PCA and generate a new 
set of models which it then distributed to all active computing 
resources. 

The authors, feeling the need to help preserve our own small 
region of one realization of the Universe, purchased carbon off- 
sets to counteract the carbon emissions resulting from this large 
investment of computing time. 
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Table 3. The set of datasets used as constraints on our model, together with a reference to where the dataset is shown in this paper 
and the value of the weight, Wi, assigned to each constraint. 



Constraint 



Reference Weight (wi) 



f4.2 



A 

10 

1 

11 

H 

14 

1 

T5 



Star formation history 

bj-band z = luminosity function F'l 

K-band 2 = luminosity function Fig, 

Morphologically segregated 2 = luminosity function Fig, 

60/im 2 = luminosity function Fig, 

Evolving K-band luminosity function Fig, 

2 = 3 UV luminosity function Fig, 

2 = 5 UV luminosity functio Fig, 

2 = 6 UV luminosity functio Fig^ 
TuUy-Fisher relation j 4.S 

Gas-phase metallicities Fig^ 20 
Colour distributions j 4.4 

Half-light radius distributions Fig. |l8 

Disk scale length distributions Fig_ 
Supermassive black hole mass distributions 

Stellar metallicities Fig, 

Gas-to-light ratios Fig^ 
Clustering f 4.Sr 

Local Group luminosity function Fig 

Local Group satellite galaxy sizes Fig 

Local Group satellite galaxy metallicities Fig 



TO 



U.9 



25 
1^ 

27 



1.00 
2.00 
2.00 
1.00 
1.00 
1.00 
1.00 
0.75 
0.75 
2.00 
1.00 
2.00 
1.50 
2.00 
1.00 
1.00 
1.00 
1.50 
1.00 
1.00 
1.00 



space in which the best model is likely to be found. We then 
adjust parameters manually to find the final model (utiliz- 
ing our intuition of how the model will respond to changes 
in parameters). 

Interesting constraints and correlations can be seen in 
Figure [T] For example, the combination ahot-Vhot,disk is 
quite well constrained and somewhat anti-correlated (such 
that an increase in Ohot can be played off against a decrease 
in Vhot.disk- It is immediately clear, for example, that no good 
model can be found with Acxpci.disk >l-5 while Acxpci, bulge is 
much less well constrained, but must be larger then about 
1.5. 

The principal component vectors from the final set of 
36,017 models are shown in Table|4] We note here that these 
vectors are quite different from those found by [Bower et al.| 
(20101. This is not too surprising as our implementation of 
Galform is quite different from theirs and we constrain our 
model to a much broader collection of datasets. We will ex- 
amine the PGA vectors in greater detail in a future paper, 
and so restrict ourselves to a brief discussion here. Taking 
the first PGA vector for example, we see that it is dominated 
by 2cut, a* and ahot- These parameters all have strong effects 
on the faint end of luminosity functions. Luminosity func- 
tions are abundant in our set of constraints and have been 
well measured. As such, they provide some of the strongest 
constraints on the model. It can be seen that an increase 
in ahot, which will fiatten the slope of the faint of a lu- 
minosity function, has a similar effect as a decrease in a*, 
which will preferentially reduce rates of star formation in 
low-mass galaxies and so also flatten the faint end slope. 
The second PGA component shows a strong but opposite 
dependence on Sib and Acxpoi, burst- Increasing Sib results in 
more fuel for galaxy formation, while increasing Aexpei, burst 
causes material to be lost by being expelled from halos. As 
we continue to further PGA vectors the parameter combi- 
nations they represent become more complicated and diffi- 
cult to interpret — the advantage of our methodology is that 



these complex interactions can be taken into account when 
exploring the model parameter space. 

The differences between our results and those of 



Bower et al. (20101 are interesting in their own right. For 



example, Bower et al. (20101 found two "islands" of good fit 
in the supernovae feedback parameter space (Virot.disk and 
'^iiot, burst): a strong feedback island (corresponding approx- 
imately to what we find in this work) and a weak feedback 
island (which we do not find). The weak feedback island is 
ruled out in the present work as, while a good fit to the 
galaxy luminosity function can be found in it (as demon- 
strated by [Bower et al. 2010[ ), no good fit to, for example, 
galaxy sizes can be found. 



4 RESULTS 

In this section we will begin by identifying the best-fit model 
and will then show results from that model compared to the 
observational data that was used to constrain the model pa- 
rameters. With the exception of results shown in §4.12[ all 
of the data shown in this section were used to constrain the 
model and, as such, the results do not represent predictions 



of the model. (In S 4.12.1 we examine the distribution of gas 
between different phases as a function of halo mass, while 
in |4.12.2[ we explore the fraction of stellar mass in the intr- 
acluster light component of halos. The data shown in these 
comparisons were not used as constraints when searching 
for the best-fit model.) The overall best-fit model (i.e. that 
which best describes the union of all datasets) is shown by 
blue lines. Additionally, we show as magenta lines the best- 
fit model to each individual dataset (as described in the 
figure captions) for comparison. We do not claim that the 
following represents a complete census of the observational 
data that could be used to constrain our galaxy formation 
model. Instead, we have selected data which spans a range 
of physical characteristics and redshifts that we think best 
constrains the physics of our model, while remaining within 
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Figure 7. Constraints on model parameters shown as 2D slices through the multi-dimensional parameter space. In each panel, the 
colour scale indicates the value of shown by the bar above the panel, with the yellow star indicating the best-fit model. Each point 

in the plane is coloured to correspond to the minimum value of found when projecting over all other dimensions of the parameter 
space. Contours illustrate the relative number of model evaluations at each point in the plane — from lightest to darkest line colour 
they correspond to 10, 30, 100 and 300 evaluations per grid cell. Most evaluations are carried out when the best model fits are found, 
indicating that our method is efficient in concentrating resources where good models are most likely to be found. To each side of the 
plane, the distribution of is projected over one of the remaining dimensions to show constraints on the indicated parameter. Top left 
panel: Shows the main parameters of the SNe feedback model, Vhot.disk and cthot- Top right panel: Shows critical parameters controlling 
the cooling and AGN feedback models, circhoat a^nd ctcool- Lower left panel: Show parameters of the adiabatic contraction model which 
have important consequences for the sizes of galaxies, Aac and uiac- Lower right panel: Shows parameters of the SNe feedback model 
that control the amount of material expelled from halos, Aoxpcl.disk and A^xpcl, burst- 



the limited (although substantial) computational resources 
at our disposal. 

In addition to these best-fit models, we will, where pos- 
sible, compare our current results with those from the previ- 
ous implementation of Galform described by [Bower et al.| 
(2006 1. Results from the Bower et al. (2006 1 model are shown 



by green lines in each figure. We have not included figures 
for every constraint used in this work — specifically, in many 
cases we show examples of the constraints only for a limited 
number of magnitude or redshift ranges. However, all of the 
constraints used are listed in Table [S] and are discussed in 
the text. 
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4.1 Best Fit Model 

The resulting set of best-fit parameters are listed in Table [S] 
We will not investigate the details of these results here, leav- 
ing an exploration of which data constrain which parameters 
and the possibility of alternative, yet acceptable, parameter 
sets to a future work. The best fit model turns out to be 
a reasonably good match to local luminosity data, galaxy 
colours, metallicities, gas content, supermassive black hole 
masses and constraints on the epoch of reionization, but to 
perform less well in matching galaxy sizes, clustering and 
the TuUy-Fisher relation. In addition, luminosity functions 
become increasingly more discrepant with the data as we 
move to higher redshifts. In the remainder of this section we 
will briefly discuss some important aspects of the best fit 
parameter set. 

The cosmological parameters are all close to the Wilkin- 
son Microwave Anisotropy Probe (WMAP) five-year expec- 
tations (by construction). The parameters of the gas cooling 
model are all quite reasonable: the three parameters arohcat 
and Qcooi are all of order unity as expected, aromovc is some- 
what smaller but still plausible, while the core radius acorc 
is around 22% of the virial radius. The parameters of the 
adiabatic contraction model differ from those proposed by 



Gnedin et al. (2004 1 but are within the range of values found 
by Gustafsson et al. (2006 1 when fitting the profiles of dark 
matter halos in simulations including galaxy formation with 
feedback. The disk stability parameter, td.gas is close to, al- 
beit lower than, the value of 0.9 suggested by the theoretical 



work of Christodoulou et al. (19951. The stripping parame- 
ter, Estrip, is of order unity as expected. 

The star formation parameters are reasonable, implying 
a low efficiency of star formation. The feedback parameters, 
Vhot, disk I burst are much lower than the value of 485 km/s 
required by Bower et al. (2006 1 and significantly closer to 



the value of 200 km/s adopted by [Cole et al. (2000[ ). This 
is desirable as values around 200 km/s already stretch the 
SNe energy budget. We also note that the value of ahot is 



lower than that required by Bower et al. (20061 and closer 
to the "natural" value of 2, which would imply an efficiency 
of supernovae energy coupling into feedback that was in- 
dependent of galaxy properties. The expulsion parameters. 



A, 



expel, disk|burst 5 



are close to unity as expected. 



The parameters of the merging model imply that mass 
ratios of 1:10 or greater are required for a major merger, a 
little low, but within the range of plausibility, while only 1:5 
or greater mergers trigger a burst. Minor mergers in which 
the primary galaxy has at least 34% gas by mass and at 
least 34% of its mass in a disk can also lead to bursts. 

Finally, the black hole growth parameters are quite rea- 
sonable: black holes radiate at about 9% of the Eddington 
luminosity, 5% of cooling gas reaches the black hole during 
radio mode feedback and around 0.5% of gas in a merging 
event is driven into the black hole. 

Overall, the parameters of the best fit model seem rea- 
sonable on physical grounds. Given the large dimensionality 
of the parameter space, the complexity of the model and 
the various assumptions used in modelling complex physical 
processes we would not consider these values to be either 



24 The [Bower et al. (2006l l model used a single value of for 
both gaseous and stellar disks. 



Table 5. Parameters of the best fit model used in this work and 
of the [Bower et al. (2006[ l model. Note that the best-fit model 
listed here is one that includes self-consistent reionization and 
evolution of the IGM (see j ]2.10| l and which has been adjusted 
to also produce a reasonable reionization history (see j4.11| . It 
therefore does not correspond to the location of the best-fit model 
indicated in Fig. [7| Where appropriate, references are given to 
the article, or section of this work, in which the parameter is 
described. 



Value 



Parameter This Work 



BowerOe 



Reference 



no 

Ao 

fib 

ho 

0-8 



0.284 

0.716 

0.04724 

0.691 

0.807 

0.933 



Cosmological 
0.250 
0.750 
0.04500 
0.730 
0.900 
1.000 



Gas Cooling Model 
Qroheat 2.32 1.260 

Ocooi 0.550 0.580 

Oremove 0.102 N/A 

acore 0.163 0.100 



Aa 



Adiabatic Contraction 
0.742 1.000 
0.920 1.000 



i2.7 



Star Formation 



iT7 



e* 


0.0152 


0.0029 


Cole et al. (2000k 


a* 


-3.28 


-1.50 


Cole et al. {'}l5m 



^djgas 



0.743 



Disk Stabilit' 
0.8O 



t ty 

4l 



Hiot.disk 
Vhot, burst 

"hot 
.^cxpcljdisk 



Supernovae Feedback 
358.0 km/s 485.0 km/s 



i 2.4.1 



f 2.12 



A, 



expel, burst 



328.0 km/s 
3.36 
0.785 
7.36 



485.0 km/s 
3.20 
N/A 
N/A 



f 2.12 



Ram Pressure Stripping 
0.335 N/A 



i2.12 



/cUip 
/burst 
/gas, burst 

P/PhVLTSt 



0.0214 
0.335 
0.331 
0.672 



Merging 
0.3000 
0.100 
0.100 

N/A 



! 2.9.1 



Black Hole Growth 
e, 0.0134 0.0398 

t;. 0.0163 N/A 

i^SMBH 0.0125 0.00500 



Cole et al. (2000 


Cole et 


al. 


(2000 


! 


i>.3 

2.5 





f 2.13 



.Malbon ct al. 



(2007J 



precise or accurate (which is why we do not quote error 
bars here), but to merely represent the most plausible val- 
ues within the context of the Galform semi-analytic model 
of galaxy formation. 

In addition to this overall best-fit model, we show in 
Table [6] the parameters which produced the best-fit to sub- 
sets of the data (as indicated). We caution that these mod- 
els were selected from runs without self-consistent reioniza- 
tion and also with relatively few realizations of merger trees, 
making them noisy. This means that, after re-running these 
models with many more merger tree realizations it is pos- 
sible that they will not be such good fits to the data. We 
do, in fact, find such cases as we will highlight below. Nev- 
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Table 6. Parameters of the overall best fit model compared to those of models which best fit individual datasets (as indicated by column 
labels). Parameters which play a key role (as discussed in the relevant subsections of ij4]l in helping to obtain a good fit to each dataset 
are shown in bold type. 



Parameter 


Overall 


Star Formation Rate 


n 


B 
o 

CO 


K20 LFs 


Morphological LFs 


Hi 
> 

II 


CO 
CO 

II 


Colours 


Tully-Fisher 


A 

Ao 


0.716 


0,723 


0.723 


0.717 


0.721 


0.720 


0.717 


0.721 


0.723 


0.722 


o 


0.04724 


0.0445 


U.U400 


0.04/9 


0.0471 


0.0491 


0.0452 


0.0477 


0.0482 


0.0441 


L 

no 


0.691 


0.677 


0.711 


0.714 


0.707 


0,726 


0.700 


0.703 


0.689 


0.724 


OS 


0.807 


0.799 


0.786 


0.805 


0.785 


0,788 


0.779 


0.765 


0.808 


0.783 


ris 


0.933 


0,955 


0.957 


0.939 


0.952 


0,960 


0.947 


0.946 


0.951 


0.959 


^reheat 


2.32 


2,68 


1.98 


1.47 


1.76 


2.34 


2.38 


2. 16 


2.26 


1.91 




0.550 


0,571 


2.31 


1. 12 


1.50 


2.67 


2.10 


2.81 


0.588 


1.06 


(^remove 


0.102 


0,842 


0.692 


0.133 


0.0986 


0.547 


U.OdU / 


0.228 


0.162 


0.125 




0.163 


0,128 


0. 168 


0.155 


0.0695 


0,187 


0.0515 


0. 142 


0.109 


0.216 


A 


0.742 


0,920 


0.819 


0.764 


0.780 


0,770 


0.804 


0.746 


0.880 


0.860 


Wac 


0.920 


0,792 


0.954 


0.941 


0.957 


0,989 


0.968 


0.972 


0.868 


0.817 


e* 


0.0152 


0,0695 


0.0453 


0.0375 


0.00520 


0,0281 


U.UzZo 


0.300 


0.0153 


0.0427 




-3.28 


-2,11 


-Z. 1 o 


-z.uo 


-2.05 


-3,43 


-z.yo 


-o.oo 


-0.392 


-1.69 


^d,gas 


0.743 


0,734 


0.743 


0.726 


0.829 


0,774 


0.804 


0.957 


0.808 


0.812 


1/ 

^hot,disk 


358.0 


425,0 






491.0 


411,0 


OUD.U 


0440. U 


459.0 


389.0 


^^ot, burst 


328.0 


130,0 


470,0 


413.0 


539.0 


498,0 


544.0 


488.0 


242.0 


370.0 


^hot 


3.36 


2,61 


2,81 


2.73 


3.57 


3,50 


3.53 


3.15 


2.95 


2.58 


•^expel,disk 


0.785 


0,738 


0,252 


0.920 


0.273 


0,477 


0.571 


0.266 


0.607 


0.551 


•^expel, burst 


7.36 


6,49 


7,90 


5.39 


5.23 


7,46 


6.62 


9.23 


6.55 


2.13 


^strip 


0.335 


0,951 


0,696 


0.248 


0.207 


0,0997 


0.739 


0.355 


0.145 


0.101 


/cllip 


0.0214 


0,184 


0,0946 


0.0658 


0.0250 


0,327 


0.0246 


0.118 


0.0250 


0.308 


/burst 


0.335 


0,260 


0,310 


0.477 


0.297 


0,286 


0.281 


0.451 


0.183 


0.263 


/gas, burst 


0.331 


0,209 


0,0817 


0.0553 


0.164 


0,349 


0.236 


0.225 


0.452 


0.160 




0.672 


0,538 


0,889 


0.890 


0.517 


0,367 


1.00 


0.215 


0.928 


0.409 




0.0134 


0,0437 


0,00596 


0.0363 


0.00877 


0,0542 


0.00423 


0.0407 


0.0130 


0.0857 




0.0163 


0.0596 


0,00476 


0.00711 


0.00188 


0,0137 


0.00728 


0.0307 


0.00788 


0.0893 


F. 


0.0125 


0,00818 


0,00289 


0.00628 


0.0206 


0,0233 


0.00190 


0.0256 


0.00271 


0.0164 


Vcut 


N/A 


17.0 


36,5 


28.4 


38.7 


43,9 


32.9 


27.7 


34.5 


12.7 


Zcut 


N/A 


10,1 


11,7 


10.9 


12.7 


12,4 


10.8 


11.9 


12.8 


10.2 



ertheless, we will refer to this table in the remainder of this 
section when exploring the ability of our model to match 
each dataset. We also point out that there is no guarantee 
that any of these models that provide a good match to an in- 
dividual dataset are good matches overall — for example, the 
model which best matches galaxy sizes may produce entirely 
unacceptable z = luminosity functions. 



4.2 Star Formation History 

Figure [8] shows the star formation rate per unit volume as a 
function of redshift, with symbols indicating observational 
estimates and lines showing results from our model. Dotted 
and dashed lines show quiescent star formation in disks and 
bursts of star formation respectively, while solid lines indi- 
cate the sum of these two. The quiescent mode dominates at 
all redshifts, although we note that at high redshifts model 
disks are typically unstable and undergo frequent instabil- 
ity events. These galaxies may therefore not look like typical 
low redshift disk galaxies. The best fit model is in excellent 
agreement with the star formation rate data from z — 1 to 



z = 8, reproducing the sharp decline in star formation below 
z = 2 while maintaining a relatively high star formation rate 
out to the highest redshifts. Our model lies below the data 
at z <.l despite being a good match to the bj-band lumi- 



nosity function (see f4.3l. This suggests some inconsistency 
in the data analysis, perhaps related to the choice of IMF 
or the calibration of star formation rate indicators. Indeed, 
the model which best fits this particular dataset (shown as 
magenta lines in Fig. [sl does so by virtue of having a large 
value of e* (see TableTGj this increases star formation rates 
overall) and a small value of acooi (which alters the crit- 
ical mass scale for AGN feedback and thereby delays the 
truncation of star formation at low redshifts). While these 
changes result in a better fit to the star formation rate, they 
produce very unacceptable fits to the luminosity functions 
(which have too many bright galaxies) and galaxies which 
are far too depleted of gas. 



The Bower et al. (2006 1 model has a much lower star 
formation rate density than our best-fit model at z > 0.5, 
although it shows a comparable amount of star formation 



in bursts. (The Bower et al. (2006 1 model still manages to 
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Table 6. (cont.) Parameters of the overall best fit model compared to those of models which best fit individual datasets (as indicated 
by column labels). Parameters which play a key role in helping to obtain a good fit to each dataset are shown in bold type. 



in Is 



Parameter 


Overall 


Tully-Fisher 


Sizes 


Metallicity 


ffl 
1.;] 
\ 

§ 


Clustering 


SMBHs 


Local Group 


Local Group 


Local Group 


Ao 


0.716 


0.722 


0.720 


0.724 


0.715 


0,714 


0.716 


0.722 


0,718 


0.722 


Ob 


0.04724 


0.0441 


0.0447 


0.0453 


0.0458 


0.0470 


0.0437 


0.0475 


0,0492 


0.0467 


ho 


0.691 


0.724 


0.698 


0.720 


0.711 


0,688 


0.699 


0.710 


0,682 


0.685 




0.807 


0.783 


0.809 


0.775 


0.788 


0.795 


0.778 


0.771 


0,769 


0.773 


ris 


0.933 


0.959 


0.961 


0.948 


0.935 


0.957 


0.945 


0.938 


0,933 


0.942 


'^reheat 


2.32 


1.91 


2.33 


2.37 


2.52 


0,922 


2.96 


2.43 


2,62 


1.81 


*^cool 


0.550 


1.06 


0.0955 


2.11 


1.48 


0,855 


1.28 


2.30 


2,25 


1.49 


O^remove 


0.102 


0.125 


0.917 


0.334 


0.146 


0,466 


0.848 


0.0814 


0,0825 


0.0508 


Q-core 


0.163 


0.216 


0.0905 


0.0772 


0.0281 


0,105 


0.222 


0.127 


0,0940 


0.0210 




0.742 


0.860 


0.964 


0.765 


0.766 


0,795 


0.876 


0.741 


0,736 


0.737 




0.920 


0.817 


0.809 


0.945 


0.989 


0,871 


0.919 


0.908 


0,928 


0.985 


e* 


0.0152 


0.0427 


0.00735 


0.00272 


0.00329 


0,0295 


0.0420 


0.00751 


0,0322 


0.0175 


a* 


-3.28 


-1.69 


-2.83 


-3.60 


-3.07 


-2,65 


-2.51 


-3.32 


-2,65 


-1.52 


^d,gas 


0.743 


0.812 


0.716 


0.774 


0.773 


0,736 


0.743 


0.957 


0,784 


0.800 


^hot,disk 


358.0 


389.0 


341.0 


497.0 


449.0 


353,0 


393.0 


374.0 


452,0 


543.0 


^^ot, burst 


328.0 


370.0 


125.0 


498.0 


496.0 


341,0 


271.0 


507.0 


533,0 


467.0 




3.36 


2.58 


3.12 


3.32 


3.53 


2,37 


3.18 


3.25 


3,14 


2.48 


-^expeljdisk 


0.785 


0.551 


0.412 


0.283 


0.380 


1,06 


0.646 


0.438 


0,659 


0.622 


-^cxpcl, burst 


7.36 


2.13 


5.62 


8.97 


7.87 


7,24 


9.86 


9.60 


8,16 


6.38 


^strip 


0.335 


0.101 


0.607 


0.0184 


0.200 


0,288 


0.359 


0.0787 


0,975 


0.595 


J^ellip 


0.0214 


0.308 


0.360 


0.0925 


0.0204 


0,107 


0.203 


0.454 


0,0672 


0.0212 


fhuTst 


0.335 


0.263 


0.242 


0.348 


0.483 


0,239 


0.435 


0.379 


0,388 


0.436 


/"gas, burst 


0.331 


0.160 


0.0937 


0.171 


0.264 


0,361 


0.120 


0.410 


0,225 


0.450 




0.672 


0.409 


0.681 


0.734 


0.825 


0,500 


0.695 


0.545 


0,251 


0.718 


e. 


0.0134 


0.0857 


0.0232 


0.0266 


0.0914 


0,0201 


0.0560 


0.0419 


0,00481 


0.00823 


V' 


0.0163 


0.0893 


0.0588 


0.00928 


0.0912 


0,0216 


0.0248 


0.0139 


0,0119 


0.00538 


F, 


0.0125 


0.0164 


0.00970 


0.00807 


0.0293 


0,00352 


0.0287 


0.0133 


0,0279 


0.00585 


Vcut 


N/A 


12.7 


27.2 


26.5 


43.0 


42,9 


28.8 


45.5 


47,5 


35.5 


^cut 


N/A 


10.2 


9.31 


11.0 


12.5 


12,7 


11.0 


12,8 


12,9 


12.7 



obtain a good match to the K-band luminosity function at 
z — however by virtue of the fact that at z <.l, where 
much of the build up of stellar mass occurs, the two models 
have comparable average star formation rates, and because 
it uses a different IMF which results in a different mass-to- 
light ratio. Our best-fit model produces has 65% more mass 



in stars at z = than the Bower et al. (2006 1 model, but 
produces only 35% more K-band luminosity density, as will 
be shown in Fig. 10 mostly from faint galaxies.) Our best- 



fit model can be seen to be in significantly better agreement 
with the data than the Bower et al. (2006) model and nicely 
reproduces the sharp decline in star formation rate at low 
redshifts. 



4.3 Luminosity Functions 

Luminosity functions have traditionally represented an im- 
portant constraint for galaxy formation models. We there- 
fore include a variety of luminosity functions, spanning a 
range of redshifts in our constraints. 

Figures [9] and [lO] show local {z « 0) luminosity func- 
tions from the Two-degree Field Galaxy Redshift Survey 



(2dFGRS) ( [Norberg et al. 2002| bj band) and the Two- 
Micron All Sky Survey (2MASS) ( |Cole et al. 2001] K band) 
respectively together with model predictions. 

It is well established that the faint end slope of the lu- 
minosity function, which is flatter than would be naively 
expected from the slope of the dark matter halo mass func- 
tion, requires some type of feedback in order to be repro- 
duced in models. The supernovae feedback present in our 
model is sufficient to flatten the faint end slope of the local 
luminosity functions and bring it into good agreement with 
the data in the bj band, except perhaps at the very faintest 
magnitudes shown. The K-band shows an even flatter faint 
end slope and this is not as well reproduced by our model. 



Both our best-fit model and the Bower et al. (2006 \ 



model produce good fits to these luminosity functions (al- 
though our best fit model produces a break which is slightly 
too bright in the K-band, indicating that the galaxy colours 
are not quite right — see S 4.4 1 . This is not surprising of course 



as these were primary constraints used to find parameters 



for the [Bower et al. (2006[ ) model. The [Bower et al. (2006 1 
model does give a noticeably better match to the faint end 
of the K-band luminosity function (although it is far from 
perfect), due to the higher value of cthot that it adopts (see 
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Figure 8. The star formation rate per unit comoving volume in 
the Universe as a function of redshift. Red points show observa- 
tional estimates from a variety of sources as compiled by |Hopkins| 
|(2004| while magenta points show the star formation rate inferred 
from gamma ray bursts by [Kistler et al. (2009[ | . The solid lines 
show the total star formation rate density from our models, while 
the dotted and dashed lines show the contribution to this from 
quiescent star formation in disks and starbursts respectively. Blue 
lines show the overall best-fit model, while magenta lines indicate 
the best-fit model to this dataset and the green lines show results 
from the [Bower et al. (2006^ model. 
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Figure 9. The ^ = bj-band luminosity function from our mod- 
els: the solid lines show the luminosity function after dust extinc- 
tion is applied while the dotted lines show the statistical error on 
the model estimate. Red points indicate the observed luminosity 
function from the 2dFGRS l |Norberg et al. 2002| . Blue lines show 
the overall best-fit model, while magenta lines indicate the best-fit 
model to this dataset and the z = K-band luminosity function 
(see Fig. |10| note that the requirement that this model be a good 
match to the ^ = K-band luminosity function is the reason why 
the fit here is not as good as that of the overall best-fit model) 
and green lines show results from the [Bower et al. (2006^ model. 



Table [5|. Unfortunately, this large value of Ohot adversely af- 
fects the agreement with other datasets and so our best-fit 
model is forced to adopt a lower value. The important point 
here is that the [Bower et al. (2006| ) model was designed to 
fit just these luminosity functions, while the current model 
is being asked to simultaneously fit a much larger compi- 
lation of datasets. This point is further illustrated by the 
magenta lines in Figs . [Oj and [ 1 0| which show the model that 
best matches these two datasets. It achieves a flatter faint 
end slope by virtue of having quite large values of ahot and 
Ocooi. This improved match to the faint end is at the ex- 
pense of the bright end though (x^ fitting gives more weight 
to the faint end, which has more data points with smaller 
error bars). 

Figure [11| shows the 6 0^m infrared luminosity function 
from Saunders et al. (1990 1 (red points) and the correspond- 
ing model results (lines). The 60/im luminosity function con- 
strains the dust absorption and reemission in our model and 
so is complementary to the optical and near-IR luminosity 
functions discussed above. Our best fit model produces a 
very good match to the data at low luminosities — the sharp 
cut off at 10^^ Lq is artificial and due to the limited num- 
ber of merger trees which we are able to run and the scarcity 
of these galaxies (which are produced by massive bursts of 
star formation). The Bower et al. (20061 model matches 



well at high luminosities but underpredicts the number of 
faint galaxies. This is due to the higher frequency of star- 




bursts at low redshifts in the Bower et al. (2006 1 model (see 
Fig. [8|, which populate the bright end of the 60/Ltm lumi- 
nosity function. It must be kept in mind that absorption 



Figure 10. The 2 = K-band luminosity function from our 
models: the solid lines show the luminosity function after dust 
extinction is applied while the dotted lines show the statistical 
error on the model estimate. Red points indicate data from the 
2dFGRS+2MASS ( [Cole et al."200T| l. Blue lines show the overall 
best-fit model, while magenta indicate the best-fit model to this 
dataset and the z = bj-band luminosity function (see Fig. [9]| 
and green show results from the [Bower et al. (2006| ) model. 
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Figure 11. The 2 = 60/xm luminosity functions from our mod- 
els are shown by the solid lines. Red points indicate data from 
[Saunders et al. (1990[ l. Blue lines show the overall best-fit model, 
while magenta lines indicate the best-fit model to this dataset and 
the green lines show results from the [Bower et al. (2006^ model. 



and re-emission of starlight by dust is one of the most chal- 
lenging processes to model semi-analytically, and we expect 
that approximations made in this work may have significant 
effects on emission at GOfim. A more detailed study, utiliz- 
ing Grasil, will be presented in a future work. The best 
fit model to this specific dataset is a good fit to the data 
although it has somewhat too many 60^m-bright galaxies. 
This is achieved by adopting a much lower value of /gas, burst 
which lets minor mergers trigger bursts more easily. This 
increases the abundance of bursting galaxies with high star 
formation rates and fills in the bright end of the GOfim lu- 
minosity function. 

Figure [12] shows the Kg-band luminosity function from 



the K20 survey ( [Pozzetti et al. 2003| ) at z = 1.0. (The data at 
z = 0.5 and 1.5 were used as constraints also.) The model 
traces the evolution of the luminosity function quite well 
but overpredicts the abundance at all redshifts. This is in 



contrast to the Bower et al. (2006 1 model which matches 



these luminosity functions quite well. This is partly due to 
the tension between luminosity functions and the star for- 
mation rate density of Fig. |8] which would be better fit if 
the model produced an even higher star formation rate den- 
sity. This constraint forces our best-fit model to build up 



more stellar mass than the Bower et al. (2006 1 model, con- 



sequently, to overpredict the abundance of galaxies at these 
redshifts. This tension between luminosity function and star 
formation rate constraints may in part be due to the difficul- 
ties involved with estimating the latter observationally (due 
to uncertainties in the IMF calibration of star formation 



rate indicators and so on; see Hopkins & Beacom (2006 1 for 



a detailed examination of these issues). The best-fit model 
to this specific dataset successfully matches the data at all 
three redshifts. It achieves this through a combination of 
relatively high (i.e. less negative) a* and a high value of 
Ohot. Together, this combination allows for a fiatter faint- 



Figure 12. The z = 1 Ks-band luminosity function from our 
models is shown by the solid lines with dotted lines indicating 
the statistical uncertainty on the model estimates. Red points 
indicate data from [Pozzetti et al. (2003| . Blue lines show the 
overall best-fit model, while magenta lines indicate the best-fit 
model to this dataset and the green lines show results from the 
[Bower et al. (2006[ l model. 



end slope while maintaining the normalization of the bright 
end. 

In addition to these luminosity functions that include 
all galaxy types, in Fig. [13] we show the morphologically 



selected luminosity function of Devereux et al. (2009 1 over- 
laid with model results. We base morphological classification 
of model galaxies on bulge-to-total ratio (B/T) in dust- 
extinguished K-band light. We determine the mapping be- 
tween B/T and morphology by requiring that the relative 
abundance of each type in the model agrees with the data 
in the interval —23.5 < Mk — 51ogjQ/i < —23.0 but the 
morphological mix is not enforced outside this magnitude 
range. Our best-fit model reproduces the broad trends seen 
in this data — although we find that too many Sb-Sbc galax- 
ies are produced at the highest luminosities. The Bower et al.[ 
|(2006| ) gives a better match to this data overall. The best 
fit to the particular dataset (magenta lines in Fig. 13 1 has a 
relatively large value of /diip, but is not significantly better 
than our best-fit model. 

In addition to these relatively low redshift constraints 
we are particularly interested here in examining constraints 
from the highest redshifts currently observable. Therefore, 
Fig. [14] shows the luminosity function of z « 3 Lyman-break 
galaxies together with the expectation from our best fit 
model (blue line). Model galaxies are drawn from the entire 
sample of galaxies at z = 3 found in the model. The model 
significantly overpredicts the number of luminous galaxies 
even when internal dust extinction is taken into account 
(the dashed line in Fig. [14] shows the luminosity function 



without the effects of dust extinction). The Bower et al 



(2006 1 model gives a similarly bad match to this data at the 



bright-end (although is slightly better at low luminosities), 
producing too many highly luminous galaxies. The best-fit 
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Figure 13. The 2 = morphologically segregated K-band luminosity functions from our models. Points indicate the observed luminosity 
function from |Devereux et al. (2009[ l for morphological classes as indicated in each panel. Blue lines show the overall best-fit model, while 
magenta lines indicate the best-fit model to this dataset and the green lines show results from the [Bower et al. (2006| l model. 
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Figure 14. The z = 3 1700A luminosity functions from our mod- 
els are shown by the solid lines with dotted lines showing the 
statistical uncertainty on the model estimates. The dashed lines 
indicate the luminosity function when the effects of dust extinc- 
tion are neglected. Red points indicate the observed luminosity 
function from [Steidel et al.| l |1999[ circles) and |Dickinson| | |1998[ 
squares). Blue lines show the overall best-fit model, while ma- 
genta lines indicate the best-fit model to this dataset and the 
green lines show results from the [Bower et al. (2006[ l model. 



model to this specific dataset turns out to be not such a 
good fit, although it is better than either of the other mod- 
els shown. The problem here is one of noise. The models 
run for our parameter space search utilized relatively small 
numbers of merger tree realizations (to permit them to run 
in a reasonable amount of time). In this particular case, the 
model run during the parameter space search looked like a 
good match to the z m 3 Lyman-break galaxy luminosity 
function, but, when re-run with many more merger trees, it 
turned out that the apparently good fit was partly a result of 
fortuitous noise. This luminosity function is particular sen- 
sitive to such effects, as the bright end is dominated by rare 
starburst galaxies. 



W 10- 
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Figure 15. The z = 5 rest-frame 1500A luminosity function 
from our models are shown by the solid lines, with statistical 
errors indicated by the dotted lines. Red points indicate data 
from [McLure et al. (2009| ) . Blue lines show the overall best-fit 
model, while magenta lines indicate the best-fit model to this 
dataset and the green lines show results from the [Bower et al.| 
[(2006[ l model. 



Finally, at the highest redshifts for which we presently 
have statistically useful data, Fi g[15[shows rest-frame UV lu- 
minosity function at z = 5 from McLure et al. (20091. These 



highest redshift luminosity functions in principle place a 
strong constraint on the model. However, the effects of dust 
become extremely important at these short wavelengths and 
so our model predictions are less reliable. As such, these con- 
straints are less fundamental than most of the others which 
we consider. We use our more detailed dust modelling for 



the Bower et al. (2006 1 model here even though the origi- 
nal [Bower et al. (2006| ) used the simpler dust model of [Cole[ 



et al. (20001. However, as noted in S 2.14.1 



in our current 



model we ensure that high-2 galaxies which are undergoing 
near continuous instability driven bursting are observed dur- 



ing the dust phase of the burst. In the Bower et al. (2006 1 
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model shown here this is not the case — such systems are al- 
most always observed in a gas and dust free state, making 
them appear much brighter. It is clear that the treatment of 
these galaxies in terms of punctuated equilibrium of disks is 
inadequate and we will return to this issue in more detail in 
a future work. 

The best fit model again overpredicts the number 
and/or luminosities of galaxies at these redshifts. The'Bower' 



et al. (20061 model performs much worse here however — 
drastically overpredicting the number of luminous galaxies. 
The majority of this difference is due to the treatment of 
dust in bursts in our current model. Additionally, however, 
this difference simply refiects the fact that high-z constraints 
were not considered when selecting the parameters of the 



Bower et al. (2006 1 model — the improved agreement here il- 
lustrates the benefits of considering a wide range of datasets 
when constraining model parameters. The best fit model to 
these specific datasets shows a steeper decline at high lu- 
minosities and a lower normalization over all luminosities. 
Once again, the best fit here is not particularly good, for 
the same reasons that the z = 3 UV luminosity function is 
not too well fit (i.e. that the models run to search parameter 
space use relatively few merger trees, leading to significant 
noise in these luminosity functions which depend on galax- 
ies that form in rare halos). This is achieved through a com- 
bination of strong feedback (i.e. high Vhot.disk) and highly 
efficient star formation with a very strong dependence on 
galaxy circular velocity. However, this achieves only a rela- 
tively small improvement over the overall best fit model, at 
the expense of significantly worse fits to other datasets. 

4.4 Colours 

The bimodality of the galaxy colour-magnitude diagram has 
long been understood to convey important information re- 
garding the evolutionary history of different types of galaxy. 
Recently, semi-analjrtic models have paid close attention to 
this diagnostic (Croton et al. 200( power et al. 2006 1. In par- 
ticular. Font et al. (2008) found that the inclusion of detailed 



modelling of ram pressure stripping of hot gas from satellite 
galaxy halos is crucial for obtaining an accurate determina- 
tion of the colour-magnitude relation. That same model of 
ram pressure stripping is included in the present work. 

Figure [16] shows slices of constant magnitude through 



the colour magnitude diagram of Weinmann et al. (20061, 



overlaid with results from our model. The model is very 
successful in matching these data, showing that at bright 
magnitudes the red galaxy component dominates, shifting 
to a mix of red and blue galaxies at fainter magnitudes. The 
median colours of the blue and red components of the galaxy 
population are reproduced better in our current model than 
by that of |Bower et al. (2006| ), although there is clearly an 
offset in the blue cloud at faint magnitudes (model galaxies 
in the blue cloud are slightly too red). Our model repro- 
duces the colours of galaxies reasonably well, so this offset 
may be partly due to the limitations of stellar population 



synthesis models. This problem with the Bower et al. (2006 1 
model was noted by [Font et al. (200"8l ) who demonstrated 
that a combination of a higher yield of p = 0.04 in the in- 



stantaneous recycling approximation ( Bower et al. (2006 1 
assumed a yield of p = 0.02) and ram pressure stripping 
of cold gas in galaxy disks lead to a much better match 



to galaxy colours. The yield is not a free parameter in our 
model, instead it is determined from the IMF and stellar 
metal yields directly (see Fig. [6|, potentially rising as high 
as p = 0.04 after several Gyr. This is very close to the value 
adopted by Font et al. (200"8| ), and our model is able to 
produce a good match to the colours. As we will see later 



(in ^4.7 1, the Bower et al. (20061 model has more serious 



problems with galaxy metallicities which are somewhat rec- 
tified in our present model thereby helping us obtain a bet- 
ter match to the galaxy colours. The best-fit model to this 
specific dataset is a better match than our overall best-fit 
model for fainter galaxies, although it performs less well at 
brighter magnitudes. At faint magnitudes it produces a bluer 
blue-cloud which better matches that which is observed. It 
achieves this success by having a much larger value (i.e. less 
negative) of q*. This parameter controls how star formation 
rates scale with galaxy mass, with this model having less de- 
pendence than any other. This improves the match to galaxy 
colours (at the expense of steepening the faint end slope of 
the luminosity function), particularly for fainter galaxies. 



4.5 Scaling Relations 

Fitting the TuUy-Fisher relation simultaneously with the lu- 
minosity function has been a long-standing challenge for 



models of galaxy formation (see Button et al. (2008 1 and 
references therein). Figure 17 shows the TuUy-Fisher rela- 



tion from the SDSS as measured by Pizagno et al. (20071 



together with the result from out best-fit model. The model 
is in reasonable agreement with zero point, although some- 
what offset to higher velocities, and in good agreement with 
the luminosity dependence and width of the TuUy-Fisher 
relation. Our new model is a significantly better match to 



the TuUy-Fisher relation than that of Bower et al. (20061, 



which produces galaxies with rotation speeds that are sys- 
tematically too large (particularly for the brightest galax- 
ies). For example, for the most luminous galaxies shown the 
Bower et al. (2006) predicts a population of galaxies with 
circular velocities of 300-400km/s or greater — strongly ruled 
out by observations. The new model on the other hand pre- 
dicts essentially no galaxies in this velocity range. The best- 
fit model to this particular dataset is a significantly better 
match than our overall best-fit model. No single parameter 
is responsible for the improvement, but Acxpoi, burst plays an 
important role — it is much lower in the best-fit model to the 
TuUy-Fisher data. 



4.6 Sizes 

Figure [18] shows the distribution of galaxy sizes, split by 
morphological type and magnitude, from the SDSS ( jShen] 
et al. 2003). To morphologically classify model galaxies we 
utilize the bulge-to-total ratio in dust-extinguished "'^r-band 



light. From the K-band morphologically segregated lumi- 



nosity function (see f4.3l we find that E and SO galaxies 



are those with B/T> 0.714 for the best fit model. There 
is no convincing reason to expect this value to correspond 
precisely to the morphological selection used by |Shen et al.| 
(20031, but it is currently our best method to choose a di- 
vision between early and late types in our model. For sim- 
plicity, we employ the same morphological cut for all three 
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Figure 16. '^ ^g— " colour distributions for galaxies at z = 0.1 split by g-band absolute magnitude (see above each panel for magnitude 
range). Solid lines indicate the distributions from our models while the red points show data from the SDSS ( [Weinmann et al. 2006| |. 
Blue lines show the overall best-fit model, while magenta lines indicate the best-fit model to this dataset and the green lines show results 
from the [Bower et al. (2006| l model. Note that the magenta model is selected on the basis on more panels than are shown here. 
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Figure 17. Slices though the i-band TuUy-Fisher relation from the SDSS ( [Pizagno et al. 2007^ at constant absolute magnitude are 
shown by red points. Solid lines show results from our models with dotted lines indicating the statistical error on the model estimate. 
Blue lines show the overall best-fit model, while magenta lines indicate the best-fit model to this dataset and the green lines show results 
from the [Bower et al. (2006| l model. 



models plotted in Fig. |18[ Model results are overlaid as lines. 
Model galaxies are too large compared to the data, by fac- 
tors of about two, and the distribution of model galaxy sizes 
is too broad. This problem is more significant for the fainter 
galaxies. 

Figure |19| shows t he distribution of disk sizes from 
de Jong & Lacey (2000 1 with model results overlaid as lines. 



This permits a more careful comparison with the model as it 
does not require us to assign morphological types to model 
galaxies. Model disks are somewhat too large in all luminos- 
ity bins considered, and the width of the distribution of disk 
sizes is broader than that observed. 



The Bower et al. (2006 1 model produces galaxies which 



are systematically smaller than those in our current best-fit 
model at bright magnitudes, but larger at faint magnitudes. 
It also produces a narrower distribution of disk sizes. Our 
best fit model to these combined size datasets is a rather 
poor match to the distribution of disk sizes. We find that 
it is challenging to obtain realistic sizes for disks in our 
model while simultaneously matching other observational 
constraints. This problem, which may reflect inaccuracies 
in the angular momentum of cooling gas, angular momen- 
tum loss during cooling or merging, or internal processes 
which transfer angular momentum out of galaxies, will be 
addressed in greater detail in a future work. 



4.7 Gas and Metal Content 

The star formation and supernovae feedback prescriptions in 
our model can be constrained by measurements of the gas 
and metal content of galaxies. Figure [20] shows the distri- 



bution of gas-phase metallicities from the SDSS (Tremonti 



et al. 20041 compared with results from our best-fit model. 



Model galaxies are drawn from the entire population of 
galaxies at z = 0.1. Tremonti et al. (20041 select star form- 



ing galaxies — essentially those with well detected H/3, Ha 
and [Nil] A6584 lines — and also reject galaxies with a signifi- 
cant AGN component. We have not attempted to reproduce 
these observational selection criteria her but note that 
excluding galaxies with very low star formation rates makes 
negligible difference to our results. The model clearly pro- 
duces a strong trend of increasing metallicity with increasing 
luminosity, just as is observed, although the relation is some- 
what too steep, resulting in metallicities which are around 
a factor of two too low at the lowest luminosities plotted. 
This relation is driven, in the model, by supernovae feed- 
back: in low luminosity galaxies feedback is more efficient 



Both because we can not, at present, include the AGN com- 
ponent in the spectra and because it would involve constructing 
mock catalogues which is too expensive during our parameter 
space search. 
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Figure 18. Distributions of galaxy half-light radii (measured in the dust-extinguished face-on r-band light profile) at z = 0.1 segregated 
by r-band absolute magnitude and by morphological class. Solid lines show results from our models while dotted lines show the statistical 
error on the model estimates. Red points data from the SDSS jShen et al. 2003| . Blue lines show the overall best-fit model, while magenta 
lines indicate the best-fit model to this dataset and the green lines show results from the [Bower et al. (2006[l model. 




Figure 19. Distribution of disk scale lengths for galaxies at z = segregated by face-on I-band absolute magnitude. Solid lines show 
results from our models while dotted lines indicate the statical uncertainty on the model estimates. Red circles show data from |de Jong| 
I&: Lacey (2 000) with upper limits indicated by red triangles. Blue lines show the overall best-fit model, while magenta lines indicate the 
best-fit model to this dataset and the green lines show results from the [Bower et al. (2006^ model. 



at ejecting material from a galaxy making it less efficient at 
self-enriching. The trend is somewhat steeper in the model 
than is observed and therefore underpredicts the metallic- 
ity of low luminosity galaxies. The spread in metallicity at 
fixed luminosity is larger than that which is observed. The 
best fit model to the metallicity datasets presented in this 
subsection can be seen to actually be a worse fit to the gas 
phase metallicity, a consequence of tensions between fitting 
this data and stellar metallcities and gas fractions. 

Figure[2T]shows distributions of mean stellar metallicity 
in various bins of absolute B-band magnitude. Data, shown 



data are quite noisy, there is, in general, good agreement of 



by points, are taken from Zaritsky et al. (19941, while re- 
sults from our best-fit model are shown by lines. For model 
galaxies, we plot the luminosity-weighted mean metallicity 
of all stars (i.e. both disk and bulge stars). Although the 



the model with this data. The Bower et al. (2006 1 model 
fails to match the scaling of metallicity with stellar mass 
seen in these data. An increase in the yield in this model 
(from p = 0.02 to p = 0.04 as required to better match 
galaxy colours; [Font et al. 2008) would improve this situa- 
tion significantly, but some reduction in the dependence of 
SNe feedback on galaxy mass is likely still required to obtain 
the correct scaling. 

Finally, Fig.[22]shows the distribution of gas-to-light ra- 
tios from a compilation of data compared to results from our 
best-fit model. Model galaxies are selected to have bulge-to- 
total ratios in B-band light of 0.4 or less and gas fractions of 
3% or more in order to attempt to match the morphological 
selection (Sa and later types) in the observations. The results 
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Figure 21. Distributions of mean stellar metallicity at different slices of absolute magnitude. Red points show observational data 
compiled by jZaritsky et al. (1994[ |. Solid lines indicate results from our models while dotted lines show the statistical uncertainty on the 
model estimate. Blue lines show the overall best-fit model, while magenta lines indicate the best-fit model to this datasct and the green 
lines show results from the [Bower et al.~(2006| ) model. 



are somewhat sensitive to the morphological criteria used, a 
fact which must be taken into account when considering the 
comparison with the observational data. The model ratio is 
somewhat too high (too much gas per unit light), but dis- 
plays approximately the correct dispersion. The |Bower et al.] 
(2006 [ ) model gets closer to the observed mean for bright 
galaxies, but shows a dramatic downturn at low luminosi- 
ties (a result of its very strong supernovae feedback). The 
best fit model to this specific dataset is an excellent match 
to both the mean and dispersion in the gas fraction data. 
This is achieved primarily via a very low efficiency of star 
formation (allowing gas fractions to stay high) coupled with 
strongly velocity dependent feedback which helps obtain the 
measured slope in this relation. 

Overall, the [Bower et al. (2006[ ) performs much less 
well in matching metallicity and gas content properties. 
This problem can be traced to the very strong scaling of 
supernovae feedback strength with galaxy circular veloc- 



ity adopted in the Bower et al. (2006 1 model and the low 
yield. This strongly suppresses the effective yield in low mass 
galaxies, resulting in them being too metal poor, and like- 
wise strongly suppresses the gas content of those same low 
mass galaxies. These constraints are among the primary 
drivers causing our best fit model to adopt a lower value 
of aiiot- 



4.8 Clustering 

Galaxy clustering places strong constraints on the occu- 
pancy of galaxies within dark matter halos and, therefore, 
the merger rate (amongst other things). To compute the 
clustering properties of galaxies we make use of the fact that 
halo occupation distributions are naturally predicted by the 
Galform model. We therefore extract halo occupation dis- 
tributions directly from our best fit model. We then employ 
the halo model of galaxy clustering ( jCooray fc Sheth 2002[ ) 
to compute two-point correlation functions in redshift space. 
These are compared to measured redshift-space correlation 
functions from the 2dFGRS ( [Norberg et al. 2002[ ) in Fig.pS] 
There is excellent agreement between the model and 
data on large scales (where the two halo term dominates). 
On small scales, in the one halo regime, the model systemat- 
ically overestimates the correlation function. This discrep- 
ancy, which is due to the model placing too many satel- 



lite galaxies in massive halos, has been noted and discussed 
previously by Kim et al. (2009| ) . In their study, Kim et al 
|(2009) demonstrated that this problem might be resolved 
by invoking destruction of satellite galaxies by tidal forces 
and by accounting for satellite-satellite mergers (both pro- 
cesses reduce the number of satellites). The current model 
includes both of these processes and treats them in a signif- 
icantly more realistic way than did |Kim et al. (2009| ). We 
find that they are not enough to bring the model correla- 
tion function into agreement with the data on small scales 
(although they do help), in our particular model. This may 
indicate that these processes have not been modelled suffi- 
ciently accurately, or that our model simply begins with too 
many satellites. We note that the [Bower et al. (2006[ ) model 
performs similarly well on large scales and somewhat better 
on small scales (the stronger feedback in this model helps 
reduce the number of satellite galaxies of a given luminosity 
in high mass halos), although it still overpredicts the small 
scale clustering, as has been noted by Kim et al. (20091. 
The best-fit model to the clustering data alone is not very 
successful. This is again due to the difficulty of computing 
accurate correlation functions using the relatively small sets 
of merger trees that we are able to utilize for parameter 
space searches, and serves as an excellent example of the 
need to include better estimates of the model uncertainty 
(i.e. the variance in predictions from the model due to the 
limited number of merger trees utilized) when computing 
goodness of fit measures. 



4.9 Supermassive Black Holes 

The inclusion of AGN feedback in semi-analytic models of 
galaxy formation necessitates the inclusion of the supermas- 
sive black holes that are responsible for that feedback. As 
such, it is important to constrain the properties of these 
black holes to match those that are observed. Figure [24] 
shows the distribution of supermassive black hole masses in 
three slices of galaxy bulge mass. Points show observational 
data from Haring & Rix (2004) while lines show results from 



our best-fit model. The model is in excellent agreement with 



the current data. The Bower et al. (20061 model produces 
nearly identical results for the black hole masses. This is 
not surprising since, as pointed out by [Bower et al. (2010[ ), 
the F, parameter can be adjusted to achieve a good fit here 
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Figure 23. Rcdshift space two-point correlation functions of galaxies selected by their bj band absolute magnitude. Solid lines show 
results from our models while red points indicate data from the 2dFGRS ( [Norberg et al. 2002| . Model correlation functions are computed 
using the halo model of clustering | |Cooray &: Sheth 2002[ l with the input halo occupation distributions computed directly from our 
best-fit model. Blue lines show the overall best-fit model, while magenta lines indicate the best-fit model to this dataset and the green 
lines show results from the [Bower et al. (2006| l model. 
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Figure 24. The distribution of supermassive black hole mass in three slices of galaxy bulge mass. Data are taken from [Ha.ring fc Rix| 
|(2004[ | and are shown by red points. Solid lines indicate results from our models with dotted lines showing the statistical uncertainty on 
the model estimate. Blue lines show the overall best-fit model, while magenta lines indicate the best-fit model to this dataset and the 
green lines show results from the [Bower et al. (2006[ l model. 



without significantly affecting any other predictions. For this 
same reason, the best fit model to these black hole data in 



not significantly better than either the Bower et al. (2006 1 
or the overall best fit model. 



4.10 Local Group 

The recent discovery of several new satellite galaxies of the 
Milky Way has lead to their abundance and properties being 
more robustly known and therefore acting as a strong con- 
straint on models of galaxy formation and has attracted sig- 



nificant attention recently ( Bullock et al. 200C penson et al. 
I2002g,3omcrvillc 2002 Gtncdi n fc Kravt sov 200eMadau et al. 
[2008e.Ma_dau ct al. 200_8haoz et al. 200£ Bov ill & Ricotti 
200^pusha et al. 200^|VIaccio et al. 2009[ ). Our model is the 



only one of which we are aware that follows the formation of 
these galaxies within the context of a self-consistent model 
of the IGM and the global galaxy population which fits a 



broad range of experimental constraints on galaxies and the 
IGM. 

To compute the expected properties of Milky Way satel- 
lites in our model we simulate a large number of dark mat- 
ter halos with masses at z = in the range 2 x 10^^- 
3 X 10^^ h~ ^ Mq. From these, we select only those halos with 
a virial velocity in the range 125-180km/s (consistent with 
recent estimates; Dehnen et al. 20Qi K.ue et al. 2008 1 and 



which contain a central galaxy with a bulge-to-total ratio 
between 5 and 20% to approximately match the properties 
of the Milky Way. This step is potentially important, as it 
ensures that the satellite populations that we consider are 
consistent with the formation of a Milky Way-like galaxjj^ 



The merging history of a halo will affect both the properties 
of the central galaxy and the population of satellite galaxies. By 
selecting only halos whose merger history was suitable to pro- 
duce a Milky Way we ensure that we are looking only at satellite 
populations consistent with the presence of such a galaxy. 
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Figure 20. Gas-phase metallicity as a function of absolute mag- 
nitude from the SDSS | |Tremonti et al. 2004[ | is shown by the 
red points. Points show the median value, while error bars indi- 
cate the 2.5, 16, 84 and 97.5 percentiles of the distribution. Lines 
indicate results form our best-fit model. Solid lines indicate the 
median model relation, dashed lines the 16 and 84 percentiles and 
dotted lines the 2.5 and 97.5 percentiles, corresponding to the er- 
ror bars on the data. Blue lines show the overall best-fit model, 
while magenta lines indicate the best-fit model to this dataset and 
the green lines show results from the [Bower et al. (2006| l model. 
(Note that dashed and dotted lines are shown only for the best-fit 
model for clarity.) 



In practice, we find that the morphological selection has lit- 
tle effect on the satellite luminosity function. However, the 
selection of suitable halos based on virial velocity produces a 
significant reduction (by about a factor of 2) in the number 
of satellites compared to the common practice of selecting 
halos with masses of approximately 10^'^/i~^Mq. Halo se- 
lection is clearly of great importance when addressing the 
missing satellite problem. We prefer to use a selection on 
halo virial velocity here rather than a selection on galaxy 

for ex- 



stellar mass, as was used by Benson et al. (2002a 



ample, since we know that the TuUy-Fisher relation in our 



model is incorrect (see [ 4.5 1 and so selecting on galaxy mass 



would result in an incorrect sample of halo masses. 

Figure [25] shows the V-band luminosity function of 
Milky Way satellite galaxies from our best fit model com- 
pared with the latest observational estimate. Our model is 
able to produce a sufficient number of the brightest satel- 
lites in a small fraction of realizations, although the median 
lies below the observed luminosity function for the Milky 
Way. At lower luminosities our best fit model overpredicts 
the observed number of satellites by factors of up to 5. It 
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Figure 22. Gas (hydrogen) to B-band light ratios at 2 = as a 
function of B-band absolute magnitude. The solid lines show the 
mean ratio from our models while the dotted lines show the dis- 
persion around the mean. Red points show the mean ratio from a 
compilation of data from Huchtmcicr fc Richte r (1988[ l and |Sage| 
[(1993) with error bars indicating the dispersion in the distribu- 
tion. Blue lines show the overall best-fit model, while magenta 
lines indicate the best-fit model to this dataset and the green lines 
show results from the [Bower et al. (2006^ model. Model galaxies 
were selected to have bulge-to-total ratios in B-band light of 0.4 
or less and gas fractions of 3% or more in order to attempt to 
match the morphological selection (Sa and later types) in the 
observations. 



has recently been pointed out ( jBusha et al. 200?^|i''ont et al.| 
2009| that inhomogeneous reionization (namely the reion- 



ization of the Lagrangian volume of the Milky Way halo 
by Milky Way progenitors) is an important consideration 
when computing the abundance of Local Group satellites. 



In particular, Font et al. (2009 1 find a similar level of dis- 



crepancy in the luminosity function when they ignore this 
effect (as we do here) and use a similar feedback model, but 
demonstrate that consideration of inhomogeneous reioniza- 
tion can reconcile the predicted and observed abundance of 
satellites. We do not consider inhomogeneous reionization 
here, but will return to it in greater detail in a future work. 
It must be noted, however, that this may have an impact 
on the luminosity function of Local Group satellites. The 



Bower et al. (2006 1 model gives a reasonably good match 



to the data, producing slightly fewer satellites than are ob- 
served at all luminosities. The best fit model to this specific 
dataset is in good agreement with the observations down to 
My — —5, but fails to produce fainter satellites. (It also pro- 
duces very few halo/galaxy pairs which meet our criteria to 
be deemed "Milky Way-like" , resulting in poor statistics for 
this model. The models utilized during the parameter space 
search happened to produce more faint galaxies, resulting 
in them being judged a good fit — this is another example 
of where understanding the model uncertainty is of crucial 
importance.) 

Figure [26] shows the distribution of half-mass radii for 
Milky Way satellites split into four bins of V-band absolute 
magnitude (only two of the bins are shown). The data are 
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enrichment in these galaxies. However, this weaker feedback 
and low ahot also result in a steeper faint end slope for the 



Figure 25. The luminosity function of Local Group satellite 
galaxies in our models. Red points show current observational 
estimates of the luminosity function from [Koposov et al. (2008| 
including corrections for sky coverage and selection probability 
from ToUerud et al. (2008). Solid lines show the median luminosity 
functions of model satellite galaxies located in Milky Way-hosting 
halos, while dotted lines indicate the lO**" and 90"' percentiles of 
the distribution of model luminosity functions. Blue lines show 
the overall best-fit model, while magenta lines indicate the best- 
fit model to this dataset and the green lines show results from the 
[Bower et al. (2006^ model. 



sparse, but the model produces galaxies that are too small 
compared to the observed satellites by factors of around 3- 
6. The Bower et al. (2006 1 model has the opposite problem. 



producing faint satellites that are too large but doing well at 
matching the sizes of brighter satellites. The best fit model 
to the Local Group size data alone is not significantly better 
than the overall best fit model — the sizes tend to be rather 
insensitive to most parameters. 

Figure [27] shows the distribution of stellar metallicities 
for Milky Way satellites split into the same four bins of V- 
band absolute magnitude (of which only two are shown). 
Once again, the data are sparse, but the model is seen to 
predict distributions of metallicity that are too broad com- 



pared to those observed. The Bower et al. (2006) model per- 
forms poorly here, significantly underestimating the metal- 
licities of the fainter satellites. This problem can be directly 
traced to the high value of Ohot used by the [Bower et al.| 
1(2006 1 ) model which results is exceptionally strong super- 
novae feedback, and consequently very low effective yields, 
for low mass galaxies. The best fit model to the Local Group 
metallicity data alone performs much better than the |Bowir| 



et al. (2006 1 and significantly better than the overall best 



fit model in reproducing both the trend with luminosity and 
scatter at fixed luminosity. This is achieved through a combi- 
nation of relatively weakly velocity dependent feedback (i.e. 
a low value of ahot) and a weak scaling of star formation ef- 
ficiency with velocity. Together, these parameters determine 
the trend of effective yield with mass and the degree of self- 



global luminosity function compared to Bower et al. (2006 1, 
thereby giving less success in matching the data in that par- 
ticular statistic. 



4.11 IGM Evolution 



As described in j ]2.10[ our model self-consistently evolves the 
properties of the intergalactic medium along with those of 
galaxies. In this section we discuss basic properties of the 
IGM (and related quantities) from our best-fit model. 

Photoheating of the IGM begins to raise its tempera- 
ture above the adiabatic expectation at z ~ 25, reaching 
a peak temperature of approximately 15,000K when hydro- 
gen becomes fully reionized before cooling to around 2,000K 
by z = 0. Hydrogen is fully reionized by z = 8. Helium is 
singly ionized at approximately the same time. There fol- 
lows an extended period during which helium is partially 
doubly ionized, but is not fully doubly ionized until much 
later, around 2 = 4. 



Figure 28 shows the Gunn-Peterson (Gunn & Peterson 



1965 1 and electron scattering optical depths as a function of 



redshift. The Gunn-Peterson optical depth rises sharply at 
the epoch of reionization becoming optically thick at z = 8. 
The rise in Gunn-Peterson optical depth is offset from that 
seen in observations of high redshift quasars, suggesting that 
reionization of hydrogen occurs somewhat too early in our 



model, although Becker et al. (20071 have argued that this 
trend in optical depth does not necessarily coincide with 
the epoch of reionization, but is instead consistent with a 
smooth extrapolation of the Lyman-a forest from lower red- 
shifts (our model does not include the Lyman-a forest). The 
electron scattering optical depth is an excellent match to 
that inferred from WMAP observations of the cosmic mi- 
crowave background (i.e. consistent within the errors) sug- 
gesting that our model reionizes the Universe at the correct 
epoch. 

One of the key effects of the reionization of the Universe 
is to suppress the formation of galaxies in low mass dark 
matter halos. We find that the accretion temperature, Tacc, 
remains approximately constant at around 30,000K below 
z = 3, corresponding to a mass scale increasing with time. 
The filtering mass rises sharply during reionization and re- 
mains large until the present day. 

We note that the model predicts too much flux at 91 2 A 
in the photon background. We suspect that this is due to 
the fact that our IGM model is uniform. Inclusion of a non- 
uniform IGM (i.e. the Lyman-a forest) would result in a 
greater mean optical depth and would reduce the model flux. 



4.12 Additional Results 

In this section we present two additional results that were 
not used to constrain the model, and therefore represent 
predictions. 



4.12.1 Gas Phases 

While not included in our fitting procedure, it is interesting 
to examine the distribution of gas between different phases 
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Figure 26. The size distribution of Local Group satellite galaxies in our models. Red points show current observational estimates of 
the size distribution from [ToUerud et al. (2008^ . Solid lines show the size distribution of model satellite galaxies located in Milky Way- 
hosting halos with dotted lines showing the statistical uncertainty on the model estimate. Blue lines show the overall best-fit model, 
while magenta lines indicate the best-fit model to this dataset and the green lines show results from the [Bower et al. (2006| ) model. 
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Figure 27. The metallicity distribution of Local Group satellite galaxies in our models. Red points show current observational estimates 
of the metallicity distribution from the compilation of [Mateo (1998| l and from [Kirby et al. (2008| . Solid lines show the metallicity 
distribution of model satellite galaxies located in Milky Way-hosting halos with dotted lines showing the statistical uncertainty on the 
model estimate. Blue lines show the overall best-fit model, while magenta lines indicate the best-fit model to this dataset and the green 
lines show results from the [Bower et al.~(2006j model. 



as a function of dark matter halo mass. Figure [29] shows the 
fraction of baryons in hot (including reheated gas), galaxy 
(cold gas in disks plus stars in disks and spheroids) and 
ejected (lost from the halo) phases. The Bower et al. (2006 1 



model (which has no ejected material) shows a peak in 
galaxy phase fraction at Mhaio « 2 x 10^^/i~^Mq with a 
rapid decline to lower mass and asymptoting to a constant 



fraction of 5% in higher mass halos. This follows the general 



trend found in semi-analytic models (see, for example, Ben- 



son et al. 2000b I in which supernovae feedback suppresses 



galaxy formation in low mass halos, while inefficient cool- 
ing and AGN feedback does the same in the highest mass 
halos. In contrast, our best-fit model shows modest ejection 
of gas in massive halos and a corresponding suppression in 
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Figure 28. Left-hand panel: The Gunn-Peterson | |Gunn fc P eterson 1965^ optical deptli as a function of expansion factor and redshift 
in our best-fit model. Points show observational constraints from Songaila (2004 blue points) and Fan et al. (200^ green points). Right- 
hand panel: The electron scattering optical depth to the CMB as a function of redshift in our best-fit model. The blue point shows the 
WMAP5 constraint punkley et al. 2009| |. 



the hot gas fraction, ahhough the trends are quahtatively 



the same as in Bower et ah (2006 1. This is different from the 
dependence of hot gas fraction on halo mass found by |Bower| 
et aL (2008 1 — our current model produces less ejection than 



found by Bower et al. (2008 1 resulting in the hot gas frac- 



tion being too high in intermediate mass halos. In particu- 
lar, the right-hand panel of Fig. |29| shows the gas fraction 
in model halos as a function of hot gas temperature. Model 
gas fractions were computed within a radius enclosing an 
overdensity of 2500, just as were the observed data. This 
radius, and the gas fraction within it, is computed using the 
dark matter and gas density profiles described in j ]2.5.4| and 
S 2.6.3 respectively. Compared to the data (magenta points), 
the Bower et al. (2006 1 model is a very poor match, show- 
ing almost no trend with temperature. Our best fit model 
also performs poorly, and it is clear that the suppression in 
hot gas fraction does not h ave the correct dep endence on 
halo masj^ In contrast, the Bower et al. (2008 1 model pro- 
duced an excellent match to these data (as it was designed 
to do). We therefore expect that our best-fit model will not 
give a good match to the X-ray luminosity-temperature rela- 
tion, and would instead require more efficient ejection, with 
a stronger dependence on halo mass in the relevant range, 
to achieve a good fit. We reiterate that these data were not 
included as a constraint when searching parameter space 
for the best-fit model. We will return to this issue in future 
work, including these constraints directly. 



Given the hot gas profile assumed in our model and the baryon 
fraction, the largest ratio of hot gas to dark matter mass we could 
find here in massive halos is 0.10 (since the gas profile is cored, 
but the dark matter profile is not). 



4-12.2 Intrahalo Light 

Stars that are tidally stripped from model galaxies become 
part of a diffuse intrahalo component which we assumes fills 
the host halo. We can therefore predict the fraction of stars 
which are found in this intrahalo light as a function of halo 
mass and compare it to measurements of this quantity. |Zi^ 
betti et al. (2005 1 have measured this quantity for clusters, 
while McGee fc Balogh (2009[ ) have measured it for galaxy 
groups. In Fig. |30| we show their results overlaid on results 
from our model. Blue points show individual model halos, 
while the blue line shows the running median of this distri- 
bution. The magenta and red points indicate the above men- 
tioned observational determinations for groups and clusters 
respectively. Our model predicts an intrahalo light fraction 
which is a very weak function of halo mass, remaining at 
20-25% over two orders of magnitude in halo mass. At fixed 
halo mass, there is significant scatter, particularly for the 
lower mass halos. Our predictions are in agreement with 
the current observational determinations, given their rather 
large errors bars, and it is clear that in the future such mea- 
surements have the potential to provide valuable constraints 
on models of tidal stripping. 



5 EFFECTS OF PHYSICAL PROCESSES 

In the previous section we have explored the effects of vary- 
ing parameters of the model and their effect on key galaxy 
properties. We will now instead briefiy explore the effects 
of certain physical processes (those which are either new 
to this work or have not been extensively examined in the 
past) on the results of our galaxy formation model. The in- 
tent here is not to assess whether these models are "better" 
than our standard model — they all utilize less realistic phys- 
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Figure 29. Left panel: Solid lines show the median fraction of baryons in different phases as a function of halo mass, while dotted lines 
indicate the lO"' and 90**^ percentiles of the distribution. Red lines show gas in the hot phase (which includes any gas in the M^ehoated 
reservoir), blue lines gas in the galaxy phase and green lines gas which has been ejected from the halo. Thin lines indicate results from 
the [Bower et al. (2006[ l model while thick lines show results from the best fit model used in this work. Right panel: The ratio of hot gas 
mass to total halo mass as a function of halo virial temperature is shown by the solid read line. Magenta points show data from |Sun| 
|et al. (2009[ l (crosses) and [Vikhlinin et al. (2009^ (squares). Both the observed data and the model results are measured within r25oo 
(the radius enclosing an overdensity of 2500). These data were not included as constraints in our search of the model parameter space. 




Figure 30. The fraction of stars which are part of the intra- 
halo light as a function of halo mass. Blue points show individual 
model halos, while the blue line shows the running median of 
this distribution. The magenta and red points indicate the obser- 
vational determinations of [McGee fc Balogh (2009l l and |Zibetti| 
|et al. (2005] for groups and clusters respectively. 



ical models — but to examine the effects of ignoring certain 
physical processes or of making certain assumptions. This 
emphasises one of the key strengths of the semi-analytic 



approach: the ability to rapidly investigate the importance 
of different physical processes on the properties of galaxies. 
Rather than showing all model results in each case, we will 
show a small selection of model results which best demon- 
strate the effects of the updated model. 



5.1 Reionization and Photoheating 

Our standard model includes a fully self-consistent treat- 
ment of the evolution of the IGM and its back reaction on 
galaxy formation. Two key physical processes are at work 
here. The first is the suppression of baryonic infall into halos 
due to the heating of the IGM by the photoionizing back- 



ground (see S 2.10.4 1. The second is the reduction in cool- 



ing rates of gas in halos as a result of photoheating by the 
same background (see |2.6.7[ ). Here, we compare this stan- 
dard model to a model with identical parameters, but with 
these two physical processes switched off. (We retain Comp- 
ton cooling and molecular hydrogen cooling, but revert to 
collisional ionization equilibrium cooling curves since there 
is no photon background in this model.) 

Figure [31] shows some of the key effects of making these 
changes to our best-fit model. In panel "a" we show the z — 
bj-band luminosity function. The model with no baryonic 
accretion suppression or photoheating (green line) shows a 
small excess of very bright galaxies relative to the best-fit 
model (blue line) due to slightly different cooling rates in this 
model which affect the efficiency of AGN feedback. As shown 
in panel "b" of Fig. |31[ the z = 5 and z — 6 UV luminos- 
ity functions are almost identical in this variant model and 
our best-fit model. At these higher redshifts AGN feedback 
has yet to become a significant factor in galaxy evolution. A 
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Figure 31. Comparisons between our best-fit model (blue lines) and the same model without the effects of suppression of baryonic 
accretion or photoionization equilibrium co olin g (green lines). Panel a: The z = bj-band luminosity function as in Fig. [9] Panel b: The 
z = 5 I500A luminosity function as in Fig. |l5| Panel c: The mean star formation rate density in the Universe as a function of redshift 
as in Fig. |8] Panel d: The luminosity function of Local Group satellite galaxies as in Fig. |25| 



small excess of galaxies is seen in the model with no bary- 
onic accretion suppression or photoheating at the faintest 
magnitudes plotted. This is as expected — those mechanisms 
preferentially suppress the formation of very low mass galax- 
ies. 

The effects of this change in the AGN feedback can be 
seen also in panel "c", where we show the star formation 
history of the Universe. At high redshifts, the two models 
are nearly identical. However, below z « 1.5 when AGN 
feedback begins to come into play, the two models diverge 
(primarily due to differences in their quiescent star forma- 
tion rates — the rates of bursting star formation remain quite 
similar), due to the weakened AGN feedback in this variant 
model. 



Finally, in panel "d" , we show the luminosity function of 
Local Group satellites. There is little difference between this 
variant model and the best-fit model for satellites brighter 
than about Afy = —10 — photoheating and baryonic sup- 
pression play only a minor role in shaping the properties of 
these brighter satellites. At fainter magnitudes, the variant 
model predicts more satellites than the best-fit model — by 
about a factor of two. Suppression of baryonic accretion and 
photoheating are clearly then important mechanisms for de- 
termining the number of satellites in the Local Group, but 
other baryonic effects (namely SNe feedback) are clearly at 
work in reducing the number of satellites below the number 
of dark matter subhalos. 
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Figure 32. Comparisons between our best-fit model (blue lines) and the same model without a full hierarchy of substructures (green 
lines). Panel a: The 2 = bj-band luminosity function as in Fig. [9] Panel b: The K-band z = luminosity function of SO galaxies as 
in Fig. |13| Panel c: The redshift space two-point correlation function of galaxies with —18.5 < bj < —17.5 as in Fig. |23| Panel d: The 
luminosity function of Local Group satellite galaxies as in Fig. |25| 



5.2 Orbital Hierarchy 

In our standard model, the full hierarchy of substructures 
(i.e. halos within halos within halos. . . ) is followed (see j ]2.8[ ). 
This is in contrast to all previous semi-analytic treatments, 
in which only the first level of the hierarchy has been con- 
sidered (i.e. only subhalos, no sub-subhalos etc.). Figure [32| 
compares results from this variant model (green lines) with 
those from our best-fit standard model (blue lines). Panel 
"a" of this figure shows the z = bj-band luminosity func- 
tion of galaxies. Without a hierarchy of substructures we find 
that this luminosity function is unchanged over most of the 
range of luminosities shown. The exception is for the bright- 
est galaxies, which become slightly brighter when no hierar- 



chy of substructures is used. These galaxies grow primarily 
through merging, and this suggests therefore that includ- 
ing a hierarchy of substructures reduces the rate of merging 
onto these galaxies. At first sight, this seems counter intu- 
itive as galaxies should have more opportunity to merge as 
they pass through each level of the hierarchy. In fact, this is 
not the case. A subhalo may sink within the potential well 
of a halo and then be tidally stripped, releasing any sub- 
subhalos it may contain into the halo. These sub-subhalos 
(which become subhalos in their new host) are placed onto 
new orbits consistent with their orbital position and veloc- 
ity at the time at which their subhalo was disrupted. The 
merging timescale for these orbits plus the time they have 
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already spent orbiting with a subhalo can be longer than 
the merging timescale they would have received if they had 
been made subhalos as soon as they crossed the virial radius 
of the host halo. This is due in part to the relatively weak 



dependence of merging timescale on rc (E) in the Jiang et al 
(2008 1 fitting formulcPjand partly due to the fact that sub- 
subhalos are ejected onto relatively energetic orbits (since 
they effectively gain a kick in velocity as their subhalo no 
longer holds them in place). 

Panel "b" in Fig. [32] shows that most of the increase 
in luminosity when the orbital hierarchy is ignored occurs 
in the SO morphological class, which, in this model, makes 
up a significant part of the bright end of the luminosity 
function. Panel "c" shows that the inclusion of the orbital 
hierarchy makes little difference to the correlation function 
of galaxies. Mergers between galaxies remain dominated by 
subhalo-halo interactions, such that this new physics has 
little impact on the number of pairs of galaxies in massive 
halos. Finally, panel "d" shows the luminosity function of 
Local Group galaxies. Their numbers are slightly reduced 
when the orbital hierarchy is ignored, a direct consequence 
of the slightly increased merger rate. 



5.3 Tidal and Ram Pressure Stripping 

Our standard model incorporates both ram pressure and 
tidal stripping of gas and stars from galaxies and their hot 
gaseous atmospheres. We compare this standard model to 
one in which both of these stripping mechanisms have been 
switched off. In general, tidal stripping of stars will reduce 
the luminosity of satellite galaxies. Ram pressure or tidal 
stripping of gas from galaxies or their hot atmospheres will 
also reduce the luminosity of satellites and, additionally, may 
increase the luminosity of central galaxies (since the stripped 
gas is added to their supply of potential fuel). 

Figure[33]compares results from the model with no tidal 
or ram pressure stripping (green lines) with our standard, 
best-fit model (blue lines). In panel "a" we show the bj- 
band luminosity function. At the faintest magnitudes, the 
model without stripping shows an excess of galaxies relative 
to the standard model. This is due to low mass galaxies in 
groups and clusters being stripped of a significant fraction 
of their stars in the standard model. Conversely, the model 
without stripping produces fewer of the brightest galaxies 
(or, more correctly, the bright galaxies that it produces are 
not quite as luminous as in the standard model). This is a 
consequence of the fact the ram pressure stripping is able 
to remove some gas from low mass galaxies, making it avail- 
able for later accretion onto massive galaxies, allowing those 
massive galaxies to grow somewhat more luminous. In panel 
"b" we examine the colour distribution of faint galaxies. The 
model with no stripping produces a shift of galaxies to the 
blue cloud as expected — with stripping included these galax- 
ies lose their gas supply and quickly turn red. 

A further effect of stripping can be seen in panel "c" 



We note that this formula has not been well-tested in the 
regime in which we are employing it. A more detailed study of 
the merging timescales and orbits of sub-subhalos is clearly war- 
ranted. 



which shows the correlation function of faint galaxies. With- 
out stripping, this is increased on small scales since a greater 
number of galaxies in massive halos now make it in to the 
luminosity range selected. Tidal stripping of stars (and, to 
some extent, ram pressure removal of gas) reduce the lumi- 
nosities of cluster galaxies and thereby reduce the number 
of galaxy pairs on small scales in a given luminosity range, 
thereby helping to reduce small scale correlations. Finally, 
we show in panel "d" the gas to light ratio in a model with- 
out stripping. In low mass galaxies the resulting ratio is 
much higher than in our standard direct result of 

this gas no longer being removed by ram pressure forces. In 
more massive galaxies there is, instead a reduction in the gas 
to light ratio relative to the standard model arising because 
much of the gas is now locked away in smaller systems and 
so not available for incorporation into larger galaxies. 

Although not shown in Fig. [33] stripping processes have 
an effect on Local Group galaxies — in the absence of strip- 
ping there is a modest increase (by around 50%) in the num- 
ber of galaxies brighter than Mv = —10, but the total num- 
ber of galaxies is mostly unchanged. Additionally, the sizes 
of Local Group satellites are larger when stripping processes 
are ignored as expected (many of the satellites lose their 
outer portions due to tidal stripping) , while metallicities are 
mostly unaffected. 

5.4 Non-instantaneous Recycling, Enrichment 
and Supernovae Feedback 

Our standard model utilizes a fully non-instantaneous model 
of recycling and chemical enrichment from stellar popula- 
tions and of feedback from SNe. We compare this model 
with one in which the instantaneous recycling approximation 
is used and in which SNe feedback occurs instantaneously 
after star formation. In this model, cooling rates are com- 
puted from the total metallicity (rather than accounting for 
the abundances of individual elements as described in j j2.6| ) 
since we cannot track individual elements in this approxima- 
tion. We adopt a yield of p = 0.04 and a recycled fraction 
of _R = 0.39 for this instantaneous recycling model. (These 
values correspond approximately the values expected for a 
single stellar population with a Chabrier IMF and an age of 
approximately 10 Gyr.) 

Figures [34] and |35] compare the results of this model 
with our best-fit standard model. In Fig.|34[ panel "a" shows 
that, at z = the bright-end of the bj-band luminosity func- 
tion is shifted brightwards in the instantaneous model. This 
is a consequence of the increased metal enrichment in this 
model which increases cooling rates (which both increases 
the amount of gas that can cool and increases the mass scale 
at which AGN feedback becomes effective) . This trend is re- 
versed at higher redshifts for the UV luminosity function 
that we consider. Here, the luminosity function is shifted 
fainter in the instantaneous model. This effect is due to in- 
creased dust extinction in the instantaneous model (which 
is able to build up metals more rapidly, particularly at high 
redshifts and so results in dustier galaxies). 

Panel "b" shows the star formation rate density as a 
function of redshift. The instantaneous model shows a lower 
star formation rate at high redshift, and a higher rate at low 
redshift compared to our standard model. At high redshift 
this can be seen to be due almost entirely to a change in the 
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Figure 33. Comparisons between our best-fit model (blue lines) and the same model without the effects of tidal or ram pressure stripping 
of gas and stars from galaxies and their hot atmospheres (green lines). Panel a: The z = bj-band luminosity function as in Fig. |9] 
Panel b: The O ig— 0.1^ colour distribution for galaxies at 2 = 0.1 with —17 < Mo.ig < —16 as in Fig. |16| Panel c: The redshift space 
two-point correlation function of galaxies with —18.5 < bj < —17.5 as in Fig. |23[ Panel d: Gas (hydrogen) to B-band light ratios at 
2 = as a function of B-band absolute magnitude as in Fig. |22| 



rate of bursty star formation. The cause of this is rather sub- 
tle: in the non-instantaneous model gas is rapidly locked up 
into stars at high redshifts and is only slowly returned to the 
ISM of galaxies. This, coupled with somewhat reduced feed- 
back in the non-instantaneous model (since it takes some 
time for the SNe to occur after star formation happens) 
makes disks more massive and therefore more prone to insta- 



bilities (see f 2.4.1 1. The non-instantaneous model has more 



instability triggered bursts of star formation at high redshift 
and there is more gas availble to burst in those events. At 
low redshifts differences in metal enrichment in hot gas in 



the instantaneous model results in slightly less efficient AGN 
feedback and, therefore, a higher star formation rate. 

Instantaneous enrichment has a big effect on galaxy 
colours as indicated in panels "c" and "d" of Fig. |34] At 
faint magnitudes we find a somewhat better fit to the data 
in the instantaneous model (the blue and red peaks are 
more widely separated and the red peak is less populated). 
However, at bright magnitudes the instantaneous model pro- 
duces too many blue galaxies and too few red ones, resulting 
in significant disagreement with the data. 

Panel "a" of Fig. [35] shows the sizes of galaxy disks. 
Remarkably, the instantaneous models shows a much better 





Figure 34. Comparisons between our best-fit model (blue lines) and the same model using an instantaneous approximation for recycling, 
chemical enrichment and SNe feedback (green lines). Panel a: The z = bj-band luminosity function as in Fig. |9] Panel b: The star 
formation rate density as a fu nctio n of redshift as in Fig. [s] Panel c: The O lg— O lf colour distribution for galaxies at 2 = 0.1 with 
-18 < Mo.ig < -17 as in Fig. [T6| Panel d: The O lg-O lr colour distribution for galaxies at z = 0.1 with -22 < Mo.ig < -21 as in 
Fig. [16] 



match to the data than our standard mode^fl This can be 
traced to a corresponding difference in the distributions of 
specific angular momenta of disks in the two models, which, 
in turn, can be traced to the different rates of instability- 
triggered bursts at high redshifts in the two models, fn the 
non-instantaneous model these happen at a high rate. As a 
result, the low angular momentum material of these disks 
is locked up into the spheroid components. Later accretion 



•^^ It is worth noting that the [Bower et al. (2006[ l model uses 
the instantaneous recycling approximation and also does better 
at matching galaxy sizes than our current best-fit model. 



then results in the formation of disks from higher angular 
momentum material, resulting in disks that are too large. 
The stochasticity of this process likewise leads to a large 
dispersion in disk specific angular momenta and, therefore, 
sizes. In the instantaneous model the rate of instability- 
triggered bursts is greatly reduced, allowing disks to retain 
their early accreted, low angular momentum material, giving 
smaller disks with less variation in size. 

Panel "b" shows an example of the distribution of 
stellar metallicities. Stars in the instantaneous model are 
enriched to higher metallicities as expected — in the non- 
instantaneous model it takes time for stars to evolve and 
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Figure 35. Comparisons between our best-fit model (blue lines) and tiie same model witli instantaneous recycling, chemical enrichment 
and SNe feedback (green lines). Panel a: The distribution of disk sizes for galaxies in the range —20 < Mj q — SlogjQ^ < —19 as in 
Fig. |19| Panel b: The distribution of stellar metallicities for galaxies in the range —20 < Mb — 51og]^Q/i < —19 as in Fig. |21[ Panel c: 
The ratio of hydrogen gas mass to B-band luminosity as in Fig. |22| Panels d & e: The gas phase metallicity as a function of absolute 
magnitude as in Fig. |20| 



produce metals, allowing less enrichment overall. Panels "c" 
and "d" show the effects on gas content and metallicity re- 
spectively. The gas content is reduced in the instantaneous 
model and is in excellent agreement with the data. This is a 
result of the late-time replenishment of the ISM in the non- 
instantaneous model by material recycled from stars. The 



instantaneous model produces lower gas phase metallicities, 
again as a result of the lack of this late-time replenishment 
which consists of relatively low metallicity material. 




Figure 36. Comparisons between our best-fit model (blue lines) and the same model without adiabatic contraction of dark matter halos 
(green lines). Panel a: The distribution of half-light radii for Local Group satellites in the magnitude range — 15 < My < —10 as in 
Fig. |26| Panel b: The TuUy-Fisher relation for galaxies in the magnitude range —21 < Ali < —20 as in Fig. |17| 



5.5 Adiabatic Contraction 

Adiabatic contraction of dark matter halos in response to the 
condensation of baryons is included in our standard model as 
described in §2.7| In Fig.[36]we compare our standard model 
with one in which this adiabatic contraction is switched off 
such that dark matter halos profiles are unchanged by the 
presence of baryons. Such a change may be expected to re- 
sult in galaxies which are somewhat larger and more slowly 
rotating. Panel "a" shows the effects on Local Group satel- 
lite galaxy sizes. A slight increase in size is seen as expected. 
For larger galaxies, we see a similar effect. Rotation speeds 
of galaxies are less affected though — panel "b" shows a slice 
through the Tully-Fisher and indicates that switching off 
adiabatic contraction has actually had little effect on this 
statistic. 



6 DISCUSSION 

We have described a substantially revised implementation 
of the Galform semi-analytic model of galaxy formation. 
This version incorporates the numerous developments in our 
understanding of galaxy formation since the last major re- 
view of the code ( |Cole et al. 2000"| . Together with changes 
to the code to implement black hole feedback (Bower et al. 



2006 Bower et al. 20081, ram-pressure stripping (Font et al 



2008|) and to track the formation of black holes ( [Malbon 
et al. 20071, we have made fundamental improvements to 



key physical processes (such as cooling, re-ionisation, galax;y 
merging and tidal stripping) and removed a number of limit- 
ing assumptions (in particular, instantaneous recycling and 
chemical enrichment are no longer assumed). In addition 
to computing the properties of galaxies, the model now 
self-consistently solves for the evolution of the intergalactic 



medium and its influence on later epochs of galaxy forma- 
tion. 

The goals of these changes have been three-fold. Firstly, 
a prime motivation has been to remove the code's explicit 
dependence on discrete halo formation events. In the older 
code, the mass-doubling events were used to reset halo 
properties and re-initalise the cooling and free fall accre- 
tion calculations. In turn, this lead to abrupt changes in 
the supply of cold gas to the central galaxy which was of- 
ten not associated with any particular merging event in 
the haloes history. The new method avoids such artificial 
dependencies and leads to smoothly varying gas accretion 
rates in haloes with smooth accretion histories, and only 
leads to abrupt changes during sufficiently important merg- 
ing events. The new scheme explicitly tracks the energetics 
of material expelled from galaxies by feedback, and also al- 
lows the angular momentum of the feedback and accreted 
material to be self-consistently propagated through the code. 
Secondly, we have aimed to enhance the range of physical 
processes treated in the code so that it incorporates the 
full range of effects that are likely to be key in determin- 
ing galaxy properties. In particular, we now include careful 
treatments of galaxy-environment interarctions (tidal and 
ram-pressure stripping), taking into account the sub-halo 
hierarchy present within each halo; we take into account the 
self-consistent re-ionisation of the IGM and the impact that 
this has on gas supply to early galaxies; and we allow for ma- 
terial to be ejected from haloes (both by star-formation and 
AGN), broadening the range of plausible feedback schemes 
included in the model. Finally, the verison of the code de- 
scribed may be driven by accurate Monte-Carlo realisations 
of halo merger trees. This allows the uncertainty in the back- 
ground cosmological parameters to be factored in to the 
model parameter constraints. 
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We have also advanced the methodology by which we 
test the model's performance by simulataneously comparing 
the model to a wide range of observational data. In addition 
to our conventional approach of primarily comparing to lo- 
cal optical and near-IR luminosity functions, we now include 
luminosity function data covering a much greater a range of 
redshift and wavelength, the star formation history of the 
universe, the distribution of galaxies in colour-space, their 
gas and metal content, the TuUy-Fisher relation and various 
observational measurements of the galaxy size distribution. 
In addition, to these galaxy properties we also use the ther- 
mal evolution of the IGM as an additional constraint. 

The drawback of introducing additional physical pro- 
cesses is that this introduces additional parameters into the 
model. However, we now beleive that we have the tools to 
efficiently explore high dimensional parameter spaces and 
thus identify strongly constrained parameter combinations, 
and the additional model freedom is much less than the sum 
of the observational constraints. We performed an exten- 
sive search of the new model's parameter space utilizing the 
"parameter pursuit" methodology of 'Bower et al. (2010 1 to 
rapidly search the high-dimensional space. 

This allowed us to find a model which is an adequate 
description of many of the data sets which were used as 
constraints. In particular, the model is a good match to local 
luminosity functions and the overall rate of star formation 
in the Universe while simultaneously producing reasonable 
distributions of galaxy colours, metallicities, gas fractions 
and supermassive black hole masses all while predicting a 
plausible reionization history. In many of the original data 
comparisons, the model gives comparable results to |Bower| 
et al. (20061. In other comparisons (particularly, colours, 
metallicities and gas fractions) it greatly improves on the 
older model. 

Additionally, most of the model parameters have shifted 
relatively little compared to the older model. Where param- 
eters have changed significantly, it is possible to identify a 
direct cause. For example, the minimum timescale on which 
feedback material can be re-accreted by a galaxy (which is 
set by Qreheat) IS shortcr for the new model. This makes good 
sense since a fraction of feedback material is now expelled 
from the system through the new expulsive feedback chan- 
nel (see ^2.12[ ). Far from indicating a lack of progress, the 
comparability of the models is a tremendous success. We 
cannot emphasise enough how much many of the internal 
algorithms of the model have been revised: the near stabil- 
ity of the end results suggests a high degree of convergence, 
and that adding additional detailing of many aspects of the 
model is not required. 

Despite this encouraging success, significant discrepan- 
cies between the model and the data remain in many ar- 
eas. In particular, the sizes of galaxies are too large in our 
model (and there is too much dispersion in galaxy sizes). 
This may refiect a break down in certain model assumptions 
(e.g. the conservation of angular momentum of gas during 
the cooling and collapse phase), or that we are still lacking 
some key physics in this part of the model model (e.g. dis- 
sipative effects during spheroid formation; |Covington et al.] 
2008). In addition to the sizes, our model continues to pro- 



from that which is observed, despite using the latest models 



of adiabatic contraction. (We note that Button et al. (20071 



duce too many satellite galaxies in high mass halos, leading 
to an overprediction of the small scale clustering amplitude 
of faint galaxies; and predicts a TuUy-Fisher relation offset 



have demonstrated the difficulty of obtaining a match to the 
TuUy-Fisher relation quite clearly, and have advocated adi- 
abatic expansion or transfer of angular momentum from gas 
to dark matter to alleviate this problem.) Additionally, at 
high redshifts the agreement with luminosity function data 
is relatively poor, but these results are highly sensitive to 
the very uncertain effects of dust on galaxy magnitudes. 

The overall aim of this work was to construct a model 
that incorporates the majority of our current understand- 
ing of galaxy formation and explore the extent to which 
such a model can reproduce a large body of observational 
data spanning a range of physical properties, mass scales 
and redshifts. This is far from being the final word on the 
progress of this model. Numerous improvements remain to 
be made — such as the inclusion of a physics-based model of 
star formation. Nevertheless, the current version has been 
demonstrated to produce good agreement with a very wide 
range of observational data. Despite the large number of ad- 
justable parameters current observational data is more than 
sufficient to constrain this model — the good agreement with 
that data should be seen as a confirmation of current galaxy 
formation theory. 

We have not attempted, in this work, to explore in detail 
which physical processes are responsible for which observed 
phenomena. That, and an investigation of which data pro- 
vided constraints on which parameters, will be the subject of 
a future work. The parameter space searching methodology 
described in this paper is quite efficient and successful, but 
is presently limited by two factors. The first is the available 
computing time and speed of model calculations which lim- 
its how fine-grained any parameter space search can be. Fur- 
ther optimization of our galaxy formation code coupled with 
more and faster computers will alleviate this problem, but 
it will remain a limitation for the near future. The second 
limitation is our ignorance about how best to combine con- 
straints from different datasets. Some of the observational 
data that we would like to use is undoubtedly affected by 
poorly understood systematic errors. As a result it is unclear 
how a precidence should be assigned to each dataset. For 
example, given the robustness of the measurements, are we 
more interested in the class of models that accurately match 
the z — 5 luminosity, or those that perform better in clus- 
tering measurements? Ideally the model would match both 
equally well, but underlying systematic errors may make this 
impossible. Furthermore, to utilize the obervational data in 
a statistically correct way we often require more informa- 
tion (e.g. the full covariance matrix rather than just errors 
on each data point) than is available. 

The most formidable challenge, however, is to better 
understand the uncertainty in each model prediction. This 
is a combination of the variance introduced by the limited 
number of dark matter halo merger trees that we are able 
to simulate and the accuracy of the approximations made in 
computing a given property in the model. The first of these 
is relatively straightforward to estimate (for example, via 
a bootstrap resampling approach), but the second is much 
more difficult. For example, we are quite sure that calcu- 
lations of dust extinction in rapidly evolving high redshift 
galaxies are very uncertain, while calculations of galaxy stel- 
lar masses at z = are much more robust. The difficulty 
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arises in assigning a numerical "weight" to the model pre- 
dictions for these different constraints. Beyond simply mak- 
ing an educated guess, one might envisage comparing pre- 
dictions of dust extinction from our model with a matched 
sample of simulated high redshift galaxies in which the com- 
plicated dynamics geometry and radiative transfer could be 
treated more accurately. The variance between the semi- 
analytic and numerical simulation results would then give a 
quantitative estimate of the model uncertainty. The prob- 
lem with such an approach is that creating such a matched 
sample is extremely difficult and time consuming. 

In addition to these uncertainties, we should really in- 
clude uncertainties arising from non-galaxy formation as- 
pects of the calculation. Good examples of these include 
the IMF (which we are not explicitly trying to predict in 
our work, but which is uncertain and makes a significant 
difference to many of our results) and the spectra of stel- 
lar populations which have significant uncertainties in some 
regimes. Understanding these various model uncertainties 
is extremely challenging, but is crucial if serious parameter 
space searching in semi-analytic models is to take place. 

However, even in the absense of a well synthesised ap- 
proach, it is clear from the data sets we have considered 
that certain key problems remain to be tackled in order 
to produce a model of galaxy formation consistent with a 
broad range of observed data. Firstly, the sizes of model 
galaxies are too large, suggesting a lack of understanding 
of the physics of angular momentum in galaxies (see j ]2.7[ ). 
It is known that the simple energy-conserving model for 



merger remnant sizes proposed by Cole et al. (20001 sys- 
tematically overpredicts the sizes of spheroids and results in 



too much scatter in their sizes ( Covington et al. 2008 1, but it 
remains unclear how much this will affect the sizes of diskj^ 
and, furthermore, many spheroids in our model are formed 
through disk-instabilities rather than mergers — there is, as 
yet, no good systematic study of how to accurately deter- 
mine the sizes of such instability-formed spheroids. The disk- 
instability process itself has signficant consequences for the 
angular momentum content of disks and, as such, a careful 
examination of this process is called for. Secondly, despite 
the inclusion of tidal stripping and satellite-satellite merg- 
ing, the number of satellite galaxies in high mass halos seems 
to remain too high, as evidenced by the clustering of galax- 



ies (see S4.8I. Thirdly, the clear tension between luminosity 



function constraints and those from the inferred star forma- 
tion rate density must be reconciled. 

The model described in this work will provide the basis 
for further improvements to our modeling of galaxy forma- 
tion. In the near future we intend to return to the follow- 
ing outstanding issues and examine their importance for the 
constraints and results presented here in greater detail: 

• when, exactly, do disk instabilities occur and precisely 
what effect do they have on the galaxies in which they hap- 
pen; 

• improved modeling of the sizes of galaxies and how dif- 
ferent physical processes affect these sizes; 

• the X-ray properties and hot gas fractions in halos and 



how these constrain the amount and type of feedback from 
galaxies; 

• the effects of patchy reionization on Local Group galaxy 
properties and on the galaxy population as a whole; 

• the importance of the cold mode of gas accretion and 
how this affects the build up of galaxies at high redshifts 



(c.f . [Brooks et al. 2009[ ) ; 

• improved modeling of AGN feedback utilizing recent es- 
timates of jet power, spin-up rates and the effects of merg- 
ers on black hole spin and mass (Boyle et al. 200^ |3enson fc 
Babul 20091); 



• examination of physically motivated models of star for- 
mation and SNe feedback utilizing the framework of[Strmger] 
fc Benson (2007[). 



7 CONCLUSIONS 

In this paper we have presented recent developments of the 
galaxy formation model Galform. This extends the model 
presented in |Cole et al. (2000[ ) and [Bower et al. (2006[ ) 
adding many additional physical process (such as environ- 
mental interactions and additional feedback channels), im- 
proving the treatment of other key processes (including cool- 
ing, re-ionisation and galaxy merging) and removing unnec- 
essary limiting assumptions (such the instantaneous recy- 
cling approximation). 

The new code is compared to wide range of observa- 
tional constraints from both the local and distant universe 
and across a wide range of wavelengths. We navigate through 
the high dimesional parameter space using the "projection 



pursuit" method suggested in Bower et al. (2010), identify- 
ing a model that performs well in many of the observational 
comparisons. We find it impossible to identify a model that 
matches all the available datasets well and there are inher- 
ent tensions between the datasets pointing to some remain- 
ing inadequacies in our understanding and implementation. 
In particular, the model as it stands fails to correctly ac- 
count for the observed distribution of galaxy sizes and the 
observed Tuly-Fisher relation. 

Galaxy formation is an inherently complex and highly 
non-linear process. As such, it is clear that our understand- 
ing of it remains incomplete and our ability to model it im- 
perfect. Nevertheless, huge progress has been made in both 
of these areas, and we expect that progress will continue 
at a rapid pace. The model described in this work provides 
an excellent match to many datasets and is in reasonable 
agreement with many others; it represents a solid founda- 
tion upon which to base further calculations of galaxy for- 
mation. In particular, with its parameters well constrained 
by current data it can be used to make predictions for as 
yet unprobed regimes of galaxy formation. 

The present work is clearly not the last word on the 
subjects covered herein, however. In fact, we expect to con- 
stantly revise our model in response to new constraints and 
improved understanding of the physicj^ This simply re- 
flects the current state of galaxy formation theory — it is a 



Disks feel the gravitational potential of any embedded 
spheroid, so their sizes will be somewhat reduced if the sizes of 
spheroids are systematically reduced. 



We intend to maintain a "living document" describing any 
such alterations at www.galform.org, where we will also make 
available results from the model via an online database. 
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rapidly developing field about which we are constantly gain- 
ing new insight. 
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